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ABSTRACT 

A  preprocessor  is  designed  to  extract  a  set  of  features 
that  enhance  natural  clustering  by  removing  extraneous 
inf  orrr^ation .  The  design  rerroves  time  shift  and  scale 
dependence  by  taking  advantage  of  invariant  properties  of  a 
Fourier  transform  followed  by  a  Mellin  transform.  The 
preprocessor  is  realized  using  an  FFT  and  a  Mellin 
transform  with  a  conventional  error  correction  term.  The 
error  term  proves  to  be  indeterminate,  but  the  error's 
bound  is  identified  as  the  envelope  for  ^'ellin  correction 
terms.  Properties  of  the  Mellin  transform  are  employed  to 
modify  the  signal  so  that  the  error  correcting  is  no  longer 
required.  The  resulting  algorithms  are  tested  with 
variously  scaled  inputs  for  which  closed  form  solutions  are 
known.  With  a  verified  modification  in  place,  the 
preprocessor  produces  features  that  are  invariant  to 
shifting  and  scaling,  while  retaining  enough  information  to 
classify  canonic  shapes,  A  method  of  improving  performance 
is  introduced. 
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I.  INTROEUCTION  AND  BACKGROUND 


Pattern  classification  is  the  assignrrent  of  a  physical 
object  or  event  to  one  of  several  prespecified  catagories 
and  is  the  result  of  an  incomplete  theory  of  perception. 
Although  rrany  transducers  are  availa'ole  for  converting 
light,  sound,  terrperature ,  reflected  radar  signals,  etc., 
to  electrical  signals,  the  ability  of  machines  to  perceive 
cr  to  recognize  their  environment  remains  very  limited.  In 
the  structured  world  of  communcations  engineering,  signals 
are  designed  to  he  detectable  and  dif f erent iable.  A  mucb 
more  difficult  problem  presents  itself  when  sensing  an 
environment  through  a  transducer  and  recognizing  or  even 
classifying  the  elements  of  that  environment  on  the  sensed 
characteristics  of  the  transducer's  electrical  output. 
Pattern  recognition  can  be  considered  a  complex 
communications  problem  (for  example,  attempting  to  teach  a 
machine  to  decode  signals  encoded  by  nature).  It  is 
possible  to  alter  the  transducer's  output  to  facilitate 
object  classification,  but  determining  how  to  alter  that 
output  is  not  a  simple  task.  The  main  elememts  of  a 
classification  system  is  shown  in  figure  1  [1].  The 
transducer  senses,  actively  or  passively,  a  set  of 
characteristics  belonging  to  the  object.  These 
characteristics   can   never  be  a  complete  description  of  an 
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A   Classification    System 
Figure   1 
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object,  but  represent,  hopefully,  enough  information  to 
classify  the  object  as  belonging  to  one  of  a  nurrber  of 
classes.  For  Instance,  temperature  is  a  characteristic  of 
the  object  and  a  class,  but  this  feature  is  of  little  value 
unless  it  differs  in  sorre  way  from  objects  belonging  to 
other  classes,  while  the  set  must  include  characteristics 
that  are  common  among  that  class.  The  preprocessor  (or 
feature  extractor)  aims  to  reduce  the  data  by  measuring  or 
quantifying  certain  properties  that  distinguish  the  sensed 
object  as  belonging  to  one  class  and  not  to  others.  This 
can  be  done  by  discerning  key  features  or  using  the  imput 
to  generate  another  set  of  features  optimized  by  some  rule. 
The  values  of  each  of  these  features  is  then  passed  to  the 
classifier,  which  evaluates  these  features  to  assign  the 
object  to  a  class. 

With  varied  success,  machine  pattern  classification  has 
been  applied  to  a  large  range  of  problems/disciplines. 
Fields  where  it  is  particularly  comrron  Include  optical 
imagery,  acoustic  signal  processing,  radiology,  radio 
astronomy,  and  electronic  warfare,  to  name  a  few.  Work  in 
many  of  these  fields  was  reviewed  during  the  development  of 
this  thesis  and  the  results  derived  and  demonstrated  here 
may,  in  turn,  be  applicable  to  the  field  in  general.  This 
effort  has  been  directed  toward  designing  a  preprocessor  to 
produce  a  set  of  features  that  are  invariant  to  information 
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known  to  be  superfluous  to  classification,  but  that  retain 
enough  information  to  classify  an  object.  The  object  has 
been  sensed  by  a  transducer  and  has  been  represented  as  an 
eirpirically  derived,  univariate  tirre  series.  Such  a  series 
would  be  the  forrr  of  data  available  froir  range  only  radar 
return  which  is  specifically  what  the  preprocessor  is 
designed  to  handle.  Returning  to  figure  1,  the  object  has 
an  infinite  set  of  characteristics  (here  portrayed  as  an 
infinite  series  of  discrete  values).  The 
transducer/receiver  has  collected  sorne  characteristics  of 
the  target  objective  in  the  presence  of  noise.  This 
inforrration  is  relayed  as  a  set  of  discrete  sarrples  (hi) 
from  a  band  limited  signal.  The  preprocessor  is  designed  to 
determine  and  code  revelent  information  (Hj)  for  the 
classifier.  If  this  task  was  done  well,  classification 
becomes  a  trivial  problem.  On  the  other  hand,  if  the 
classifier  becomes  ideal  (capable  of  resolving  an  infinite 
number  of  characteristics  in  noise)  the  preprocessor  design 
begins  to  look  like  a  wire.  The  distinction  between  the 
preprocessor  and  the  classifier  is  arbitrary  from  an 
analytical  point  of  view.  When  designing  a  classification 
system  functionally,  a  difference  is  enforced.  The 
classifier  has  little  concern  for  how  the  features  are 
developed,  but  seeks  to  efficiently  use  those  provided  to 
guess   the   class  of  the  target  object.  The  preprocessor  is 
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problerr  dependent,  needing  to  produce   an   optirral   set   of 
features,  E j ,  from  the  sensed  data  hi. 

A.  FEATURE  EXTRACTION 

For  the  purposes  of  this  paper,  two  generic  approaches 
to  feature  extraction  are  defined.  The  first,  a 
classification  approach,  was  descrihed  ahove.  The  second,  a 
descriptive  approach,  tries  to  define  the  object  in  terms 
of  the  objects'  structural  features.  This  system  might 
recognize  a  car,  for  example,  by  breaking  up  the  visual 
picture  into  canonic  shapes,  and  comparing  this  to 
previously  specified  canonic  class  models.  The  perceived 
structure  of  the  physical  object  is  maintained  and  should 
reflect  the  structure  of  the  object  itself.  This  approach 
could  be  robust  to  temporary  changes  in  the  object  itself. 
In  the  car  example,  knowing  that  at  one  end  of  the  car  the 
trunk  can  be  opened  or  closed  allows  the  device  to  take 
this  factor  into  account.  Another  important  advantage  to 
descriptive  techniques  is  that  the  class  characteristics 
may  be  entered  or  specified  without  collecting  actual 
transducer  generated  data  to  train  the  machine. 
Unfortunately,  the  problem  of  designing  a  machine  to 
analyse  a  visual  scene  to  produce  a  structural  description 
has  proved  to  be  quite  difficult.  Object  description  from  a 
univariate   time   series  is  even  more  difficult,  and  if  the 
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radiation  sensed  ^y  the  transducer  is  not  frorr  the  visual 
spectrum,  the  task  rapidly  approaches  the  impossible.  For 
these  reasons  the  approach  taken  was  the  classification 
approach  (to  reduce  the  signal  to  a  set  of  orthogonal 
features  that  do  not  uniquely  reflect  the  structure  of  the 
object,  hut  do  retain  sufficient  information  to  classify 
the  object ) . 

This  paper  excludes  a  detailed  description  of  the 
transducer  specification.  The  problem  of  the  classifier 
Itself  is  viewed  as  one  of  partitioning  the  feature  space 
(Ej)  into  regions;  one  region  for  each  category.  Ideally, 
this  partitioning  should  be  arranged  so  that  none  of  the 
decisions  are  ever  wrong.  When  this  cannot  be  realized,  at 
least  the  probability  of  error  should  be  minimized  or  the 
average  cost  of  errors  minimized.  The  problem  is  one  within 
statistical  decision  theory,  knowledge  of  the  object 
classes  (the  transducer  and  the  classifier)  are  required  to 
design  the  preprocessor,  which  is  the  topic  of  this  thesis. 
The  preprocessor  designed  and  built  here  generates  features 
from  a  range  only  radar  video  signal.  These  features  are 
used  by  a  general  Bayesian  learning  classifier.  The 
supervised  learning  general  Bayesian  classifier  approaches 
the  problem  by  taking  a  series  of  incoming  sets  of  features 
labeled  as  to  their  class.  From  the  data,  an  a  posteriori 
density  is  computed.  Each  successive  set  of   training   data 
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is  used  to  refine  the  densities'  statistics.  When  the 
classifier  has  been  trained  on  N  classes,  the  features  are 
irodified  to  separate  the  class  volurres  in  an  optical  way 
and  to  reduce  the  nurrher  of  features  to  one  less  than  the 
number  of  possible  classes.  The  feature  vectors  of  class 
rembers  are  clijstered  about  a  siiTplei  point  and  likely- 
boundaries  are  set  up  allowing  classification  of  the  object 
as  belonging  to  one  of  the  classes,  or  of  an  unkown  class. 
An  \  siirplex  is  a  collection  of  N  points  in  (N-1)  space 
where  the  distance  between  any  two  of  the  points  is  equal 
to  the  distance  between  any  other  two.  Thus  a  three  class 
probletr  transf  orired  into  a  three  sirrplei  in  a  two 
dimensional  plane  produces  clustering  of  the  three  class's 
about  the  vertices  of  an  equilateral  triangle.  The  simplex 
coordinates  are  the  reduced  feature  vectors,  generalized 
from  the  training  data.  In  a  controlled,  simple  problem  the 
classifier  works  well,  but  when  encountering  real  problems 
one  class's  feature  space  will  intertwine  another's,  making 
it  much  more  difficult  to  obtain  separation  in  a  meaningful 
way.  The  goal  of  the  preprocessor  is  to  present  key 
features  that  determine  class  for  subsequent  optimization 
by  the  classifier. 
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B.  FOURIIR  -  MELLIN  PREPROCESSING 

Coirmon  to  all  univariate,  tirre  series  classification 
probletrs  are  several  variables  that  interfere  with  the 
recognition  process.  Assurring  discrete  data  processing  is 
used,  these  are  addressed  in  the  following  order? 
windowing,  fraining,  scaling,  sampling  rate,  quantization 
noise,  and  sufficient  information.  For  real  processing,  the 
iirput  waveform  is  not  sampled  for  all  tine.  It  is  sampled 
for  a  period  of  time.  This  windowing  of  the  data  corrupts 
the  resulting  spectrum  in  two  ways  [2].  First  it  introduces 
a  periodicity  (the  inverse  of  the  window  length)  and 
resulting  aliasing  to  the  otherwise  infinite  spectrum,  and 
further  distorts  the  spectrum  hy  a  convolution  with  the 
spectrum  of  the  window  itself.  3oth  of  these  effects  will 
color  all  of  the  data  in  the  same  way  and  so  can  he 
accounted  for  hy  deterministic  methods.'  Framing  can  he 
considered  characteristic  of  poor  synchronization,  whereby 
the  pattern  of  concern  is  not  position  stable  with  repect 
to  the  window  as  shown  in  figure  2a.  It  seems  that  even  in 
human  optical  recognition,  the  eye  tends  to  center  the 
pattern  prior  to  recognition,  unless  trained  otherwise. 
Scaling  is  that  property  whereby  the  object  field  may  vary 
in  scale  or  aspect  angle,  in  one  dimension  or  several 
dimensions  as  in  figure  2b.  Before  the  analog  signal  is 
sampled,   it   must   be  fed  through  a  sharp  low  pass  filter. 
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because  no  higher  frequency  noise  can  he  present  without 
being  folded  onto  the  valid  data.  Quantization  noise,  due 
to  the  requirement  to  round  off  each  sample  to  some 
discrete  level,  is  treated  the  same  as  round  off  error  in 
nurrerical  processing  [3].  In  all  of  these  problems 
discussed  to  this  point,  the  effect  of  this  processing  is 
to  mask  the  actual  feature  vectors,  making  the 
classification  system  less  sensitive  to  valid  pattern 
characteristics.  In  all  recognition  problems,  it  is  assumed 
that  there  is  sufficient  information  present  for  a  pattern 
to  be  detected  and  identified  by  the  system.  This  means 
that  there  is  sufficient  variability  between  the  classes, 
but  sufficient  similarity  between  those  patterns  of  a 
single  class  to  classify  each  pattern  in  terms  of  those 
classes . 

A  verifiable  goal  of  this  preprocessor  is  to  produce  a 
set  of  features  that  are  invariant  to  shifting  and  scaling 
changes.  An  approach,  figure  3,  has  been  proposed  and  used 
[4-10] .  The  preprocessor  consists  of  putting  the  sampled 
data  successively  through  two  transforms,  a  fast  Fourier 
transform  and  a  discrete  Mellin  transform.  Integral 
transforms,  the  Fourier  and  Mellin  included,  develop 
naturally  from  the  solution  of  simple  problems  in  potential 
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theory  [11-13]  .  The  Fourier  integral  transforrr  of  the 
waveform  h( t ) , 

-«» 

where  I=exp  [-2f^f t] ,  is  a  principle  analytical  tool  in  such 
diverse  fields  as  linear  systems,  optics,  probability 
theory,  quantum  physics,  and  signal  analysis  [13].  Its 
purpose  in  the  preprocessor  is  twofold,  but  relies  on  a 
single  characteristic.  The  magnitude  of  the  Fourier 
transforrr  is  invariant  to  shifting,  h(t-a). 


lH(^ij     =  /  e"'  fj(s)l 


(2) 


This  characteristic  removes  the  effects  of  framing 
inaccuracies,  and  also  permits  the  averaging  of  successive 
looks  or  pulses  of  data  to  improve  feature  resolution,  but 
removes  much  of  the  information  about  a  signal's  structure 
as  discussed  later.  A  discrete  Fourier  transform, 

Hc/;  ->  H  (^/t;    ^»  fl;  /;  ?y  -'J   ^- '       (3) 


l^(^)    =*  X  ^^'^-^  «'^ 


(4) 


yrr\zo 


has  an  equivalent  identity,  but  it  is  only  exact  for  shifts 
of  integer  values. 
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Shifts  of  other  than  Integer  values  result  in  errors  that 
depend  not  only  on  the  shift,  hut  on  qualities  of  the 
sarrpled  waveform  itself,  h(t). 

The  Mellin   transform  is  an  integral  transforir  with  the 
kernel,  K. 


K   =  t 


*-/ 


(6) 


H  c^)  «  J  lU)  t'-'cit. 


(7) 


is  the  Mellin  transform  with  respect  to  the  complex 
parameter  s=r-j2-f>m.  Several  simple  substitutions  relate 
this  to  more  common  analytical  tools.  Exponentially  warping 
t=exp[x]  changes  the  appearance  of  the  integral  to  what  is 
often  called  a  double  sided  Laplace  transform  [14], 


(8) 


The   transform   is  invariant  to  t  domain  scaling  when  taken 
with  respect  to  the  imaginary  part  alone, 


//(^; 


;  e-^^"-*^  V^ 


(9) 
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Equation  (9)  is  recognized  as  the  Fourier  intergral  of  an 
exponentially  distorted  wavefortr  h'(x).  The  rr^odulus  of  the 
expression  on  the  right  is  the  iragnitude  of  the  Fourier 
transforiT  of  the  exponentially  distorted  time  function.  The 
property  to  be  exploited  in  a  Fourier-Mellin  (FM) 
preprocessor  is  that  the  rrodulus  of  the  transform  in  s,  is 
invariant  to  t-scaling.  Given  a  time  waveform  h(t),  its 
fellin  transform  H(s)  is  given  as  equation  (7).  Scaling  the 
entire  t-domain  hy  k, 

letting  r=t/k,  and  remembering  that  s  is  imaginery, 

JA'  l^cs)j   =  /  Mcsjj  (11) 

Much  effort  and  detail  is  spent  implementing  equation  (7) 
digitally  in  Chapters  II  and  III.  The  rest  of  this  chapter 
is  devoted  to  providing  some  required  background  on  the 
discrete  versions  of  the  Fourier  transform  and  some 
properties  upon  which  the  FM  preprocessor  depends.  The 
treatment  here  is  brief,  being  mainly  a  review  of  basic  FFT 
concepts  and  as  such  may  be  skipped  without  loss  of 
content. 

A  discrete   Fourier   transform   is   not   computationally 
efficient   and   so   leads   to  imprac tically  long  processing 
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times.  The  fast  Fourier  transform  (FFT)  efficiently 
computes  the  discrete  transform  and  is  used  in  the 
preprocessors  "built  for  this  thesis.  Other  properities  of 
the  FFT  should  he  presented  hefore  getting  into  the 
detailed  design  of  the  preprocessor.  The  first  and  last  are 
concerned  with  symmetry.  If  h(m)  is  real,  as  in  the  case  of 
the  sampled  data,  then  the  frequency  spectrum  of  that  data 
is  even,  j Ee(n) 1 = ! He (-n ) , .  H(n)  has  a  real  part  and  an 
imaginary  part,  P.e(n)  and  Im(n),  while  h(m)  has  an  even 
he (m)=he (-m ) ,  and  an  odd  ho (m) =-ho (-m)  part. 

7?e  Cm;  =*  2,  vfe  C/m)  cts    \  — ;q^ j 


(12) 

The  odd  part  of  h(n)  times  the  cosine  kernel  summation,  and 
the  even  part  of  h(a)  times  the  sine  kernel  summation  are 
hoth  zero.  H(n)  can  now  be  seen  to  have  an  even  and  odd 
part.  Taking  the  magnitude  of  an  odd  function  makes  it 
even,  proving  that  the  Fourier  transform  of  a  real  series 
is  a  spectrum  whose  magnitude  is  symmetric  about  f=0.  This 
is  of  course  true  for  both  the  integral  and  discrete 
Fourier  transforms.  The  second  line  of  symrretry  is  an 
effect  of  discretizing  the  signal  and  its  spectrum  for  a 
discrete  transform  as  evinced  by  its  theoretical 
development.  Before  proceeding  though,  the  convolution 
theorem  is  required. 
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The   convolution   of   the  two  functions  h(t)  and  g(t)  is 
defined  as  the  farriliar  integral, 

The  relationship  between  convolution  and  the  Fourier 
integral  is  very  important  to  modern  analysis  and 
contributes  to  making  the  Fourier  transform  a  key  analytic 
tool.  The  convolution  theorem  states  that  if  h(t)  has  a 
Fourier  transform  F(f)  and  g(t)  has  the  Fourier  transform 
G(f),  then  h(t)*g(t)  has  the  Fourier  transform  H(f)G(f). 

Jlytj]   =  HU)<^(^)  (14) 

Proving  this,  the  Fourier  integral  is  used  directly  on  the 
convolution  integral. 


=  r**[  r  ^c^)ji(-6't-)^^i  e~^  ^  c/t 


Yc-f) 


From  the  shifting  theorem  already  presented. 


(16) 
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And    finally. 


Yc^)^JC^c-t)ie£ci:)]=  Cc/;//c/; 


(17) 


It    can   be    shown    similarly   that. 


T'L<^C^)  ^  MCf)]  =  ^  U)J<U) 


(18) 


Clearly,  convolution  in  one  domain  is  simple  multiplication 
in  the  other  domain.  Although  not  needed  for  the  pending 
development  of  the  discrete  Fourier  transform,  another 
important  relationship,  k:nown  as  the  correlation  theorem, 
can  he  appropriately  dealt  with  here.  The  correlation 
integral , 


(IS) 


has  an  operation  with  which  it  forms  a  Fourier  pair  as  did 
the  convolution-multiplication  operations.  This  theorem  can 
"be  established  as  before, 


^J^X'^k 


=  CfC^;[]  Jci')e»s  r^f^'^-^V^  y.^'  j  ^Ct)sIk  Clfl»/>)«^T  / 


(20) 
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The   term   in  brackets  is  the  complex  conjugate  of  H(f)  and 
is  denoted  by  H*(f)  in  the  final  forrr  of  the  theoren?  below. 


jlf^M^Ct^-DoL^]    -  Q(f^    ^*<^) 


(21) 


We   will   now   continue  on   to   the   discrete  Fourier 

transform   starting  with   a   continuous  waveform  hU).  The 

waveform  is  sampled  or  multiplied   by  a   string  of   delta 
functions ,  s( t) . 


00 


(22 


/m  ■-< 


where   capital  delta  is  the  sampling  interval.  The  infinite 

sum  is  not  realized  and  must  be  windowed,  in  this   exam.ple, 

by  w(t)=l   for   0  -t  5(M-1)4   and  zero  elsewhere.  So  that 

now , 

ra-i 


(23) 


/*>\»0 


Recalling  the  convolution  theorem,  the  multiplications  in 
the  time  domain  correspond  to  convolutions  in  the  frequency 
domain  with  the  following  results.  H(f)  is  convolved  with 
the  window  functions'  spectrum  and  will  have  the  apparent 
effect  of  Introducing  ripples  because  of  the  window's 
significant  sidelobes.  The  rippling  may  be  minimized  by 
choosing  a  window  function  with  small  sidelobes  at  the  cost 
of   other,   perhaps   more  acceptable,  spectral  degradation. 
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H(f)  is  convolved  with  the  sar^pling  function  and 
S(f)=I(f-n/<d  )  has  made  the  spectrum  periodic  with  respect 
to  the  interval  ?=1/A  .  The  spectrurr  coming  from  a  real 
waveform  is  first  symmetric  about  f=0  as  shown  before.  Now 
because  of  its  periodocity  the  spectrum  is  symmetric  also 
about  f=F/2  fthe  Nyquist  or  folding  frequency).  One  final 
step  remains.  The  Fourier  spectrum  is  also  taken  at 
discrete  points  1/T  apart.  The  result  in  the  time  domain  is 
the  convolution  of  the  sampled,  windowed  signal  with 
I(t-nT),  which  is  a  periodic  signal  with  T  as  its  period. 

The  FFT  algorithm  used  in  the  preprocessors  developed  in 
this  thesis  uses  a  Cooly-Tuicey ,  base  two  algorithm  [15]. 
This  is  documented  where  it  is  used  in  programs  included  in 
the  Appendices.  The  algorithm  uses  N  samples  where  N  is  two 
to  an  Integer  power.  In  the  transformed  domain,  due  to  the 
symmetry  shown  above,  there  are  N  coefficients  (only  N/2  of 
which  are  unique  as  shown  in  figure  4).  The  original 
waveform  must  be  band  limited  prior  to  sampling  to  minimize 
aliasing.  The  resulting  frequency  spectrum  should  approach 
zero  at  the  folding  frequency  where  up  to  half  of  the  power 
can  be  aliasing  noise. 

In  the  following  two  chapters,  several  different 
preprocessing  algorithms  are  developed  and  their 
performance  compared  with  canonic  inputs.  In  Chapter  IV  the 
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pro"blerr  of  ship  identification  with  range  only  radar  is 
discussed  briefly  and  one  preprocessor  is  applied  to 
several  ship  profiles.  In  the  fifth  and  final  chapter, 
conclusions  are  drawn  from  the  work  done  here  and  follow-on 
efforts  are  recommended. 
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An  FFT  Spectrum 
Figure  4 
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II.       DIGITAL   MILLIN    TRANSFORM 
BY   EXPONENTIAL    WAPPING 

In  the  first  chapter,  the  rrodulus  of  the  i^ellin 
transform  was  shown  to  be  invariant  to  scaling.  A  detailed 
eiairination  of  the  irechanics  involved  suggests  an 
Inrplementation  that  is  widely  used.  The  Mellin  transforrr  of 
a  t-doirain  function  h(t)  is  given  in  equations  (7)  and 
(9'). 

Res)  -  ^"^IC^eV^^y^  (S') 

-09 

Delta  has  heen  added,  corresponding  to  the  sample  interval 
in  the  t-domaln.  Again  it  is  noted  that  (9')  is  a  Fourier 
transform,  where  s=-J2'frm.  Solving  the  integral  for  the 
effect  of  a  t-scaling  ty  the  factor  k, 

^^^-^^Hcs)  (24) 

Clearly,  in  the  Fourier  integral,  the  scaling  factor  k  has 
become  a   shift   for  which  the  modulus  of  the  transform  is 

invariant.  The  exponential  warp  alone  has   transformed  the 

scale   factor   into   a   shift.   A  prerequisite  is  that  the 

t-domain  signal  has  no  shift.  If  there  were   a   shift,  it 

does   not   transform   to   a   simple   factor  or  shift  in  the 
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Mellin  dorrain  and  so   cannot   subsequently   "be   rerroved   by 
taking  the  magnitude  of  a  Fourier  transforrr. 

Irrplerrentation  of  a  discrete  Mellin  transform  is  as 
difficult  as  it  is  with  the  Laplace  transform  [13].  Once 
transformed,  characteristics  in  the  new  spectrum  are 
difficult  to  relate  to  the  original  signal  characteristics. 
One  hypothesis  relating  the  two  domains  is  generated  by 
comparison  to  the  Fourier  transform.  The  power  spectrum 
associated  with  the  Fourier  transform  can  be  used  to  detect 
periodicities  in  the  physical  function,  since  the  wave 
numbers  at  which  sharp  peaks  of  the  spectrum  occur  give  the 
wavelength  of  such  periodicities.  By  analogy,  the  positions 
of  the  peaks  in  the  spectrum  associated  with  the  Mellin 
transform  are  said  to  give  the  magnification  or  compression 
which  will  produce  features  in  the  physical  domain. 
Further,  this  stretching  and  compressing  is  identified  as 
periodic  in  nature  [21].  This  seems  unlikely  because  the 
Mellin  is  invariant  to  magnification/compression  and  does 
not  behave  well  (a  scale  factor  that  is  a  function  of  t  and 
k(t)  as  seen  in  (24).  The  Fourier  spectrum  models  the 
original  signal  by  a  set  of  weighted  sinusoids  of  varying 
frequencies,  and  therefore,  naturally  display  periodicity 
and   is   invariant   to  shift.  The  Mellin  is  also  a  weighted 
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surr  of  sinusoids,  but  whose  rragnitudes  are  inversely 
weighted    by    t. 

Has)  ^    \  Act)    ^^ ^±  (25) 

Values  for  h(t)  for  0  t  1  are  far  irore  important  to  the 
sum  than  those  beyond  that  point.  This  characteristic 
contributes  to  the  difficulty  of  realizing  discrete  Laplace 
and  Mellin  transforms  and  is  a  major  topic  covered  in  this 
chapter . 

The  numerical  approximation  of  the  Fourier-f^ellin  (Jt^) 
transformation  by  exponential  warping  is  functionally 
diagrammed  in  figure  5.  A  FORTRAN  program  using  the 
algorithm  described  in  this  chapter  is  included  in  Appendix 
A.  Referring  to  figure  5,  the  input  samples  are  from  a 
pulse  whose  duration  is  finite  and  less  than  that  of  the 
sample  window.  For  this  case,  no  spectral  distortion  is 
experienced  by  filling  zeros  behind  the  sampled  data.  The 
only  effect  of  the  zero  filling  is  to  add  spectral 
resolution  in  the  frequency  domain.  Once  through  the  first 
FFT  block,  and  the  magnitude  taken,  Ef  is  symmetric  about 
zero  and  at  the  folding  frequency.  Thus,  for  N  time  samples 
and  filled  zeros,  there  will  be  N/2  unique  spectral 
coefficients.  This  unique  spectrum  is  interpolated, 
resampled  as  a  warped  signal  and  fed  to  the  final  FFT 
block.  A  correction  is  added,  and  the  modulus  is  taken.  The 
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resulting  FM  features  are  invariant  to  shifting  and  scaling 
in  the  tiire  domain.  The  FJT  block  was  covered  in  sufficient 
detail  in  the  preceding  chapter.  Some  effort  will  be  spent 
in  discussing  the  warping  itself  and  the  need  for,  and  the 
development  of,  a  zero  point  correction. 

A.  ALGORITHM  DEVELOPMENT 

This  section  uses  its  own  notation  to  address  the 
requisite  exponential  sarrpling.  The  series  to  be 
transformed  is  h(f).  Its  Mellin  transform  is  E(m)  where  m 
is  the  Mellin  frequency  in  s=-j2'l^m.  The  transform  is, 

Hc-)=  f*iW/''W^  (26) 

Letting  f=Fexp[x]  as  before,  where  F  is  added  corresponding 
to  the  sample  interval 

The  need  is  to  evaluate  M  equally  spaced  samples  in  x,  at 

^^X,  :^x,  ..  .  ,  C  M-n  X  (28) 

while  the  data  consists  of  N  equally  spaced  samples  in  f, 

0,  F,  zr,    ...  ,  CN-i^  F  (29) 

Assuming  that  F  is  a  small  enough  interval  to  properly 
characterize  h(f )  ,   it   is   sometimes   advisable  [9,10]  to 
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choose  the  exponential  sampling  interval  (X)  such  that  the 
largest  intersample  spacing  in  h(rexp[xl)  is  equal  to  IcF 
where  k=l .  The  other  set  of  conditions  used  to  uniquely 
specify  the  new  samples  are  that  f=F  and  x=0  will  he  the 
lowest  lirrit  for  interpolation,  while  (N-l)F  and  (r^-DX  are 
equated  as  the  upper  limit.  The  first  requirerrent 
constrains   the   choice   of   M  by 


L  =    e^*^"'^'^~     ^c'M-z)A 


(30) 


Meeting   the   second   condition,   the   end   points   in  each 
sampled  series  are  equated  yielding. 


C/y-/^  =1  e 


c<vj-f;->f 


(31) 


Substituting  (31)  into  (30)  while  applying  the  exponential 
series  approximation  gives, 

y=k   -^  (32) 

^  N-/ 

One  more  substitution,  (32)  back  into  (31)  sets  up  the 
desired  result  where  M  is  now  specified  to  exponentially 
sample  from  f=F  to  f=(N-l)F  with  kF  being  the  largest 
interval  between  samples. 


J! 


(33) 


34 


As  N  gets  larger,  M  explodes.  If  N=ie  then  ^'=41 ,  if  N=32 
then  M=106,  if  N=64  then  M=261,  and  so  on.  This  strict 
requlreirent  can  be  corrprorrised  depending  upon  the 
application.  The  other  extreme  [18]  is  the  criterion  that 
requires  the  strallest  interval  "between  samples  to  equal  the 
interval  between  uniformly  spaced  Nyquist  samplinc^,  with 
the  intervals  increasing  exponentially  thereafter.  The 
specification  permits  analyzing  frequencies  approaching  the 
largest  which  can  be  analyzed  with  uniform  spacing.  This 
greatly  reduces  the  number  of  samples  and  decreases  the 
required  processing  tirre,  but  is  limited  in  application. 
Two  factors  mitigate  the  stringent  requirements  imposed  by 
(33)  where  k=l .  First,  using  N  unique,  uniform  samples 
yields  N/2  unique,  uniform  samples  in  the  FTT  domain.  A 
related  consideration  is  that  the  values  of  the  spectral 
components  approach  zero  at  the  folding  frequency  because 
the  original  signal  has  been  band  limited  and  some  over 
sampling  is  normally  recommended.  Secondly,  the  inverse  f 
weighting  apparent  in  equation  (25) 

H(^)i  I    ±1^   e  ^       Uf  (25') 

attaches  a  decreased  importance  to  h(f)  near  the  folding 
frequency  fn.  These  two  effects  combine  to  permit  a  much 
lower  sampling  rate  once  the  modulus  of  the  FFT  has  been 
taken.   In  spite  of  this,  the  stiffer  rule  (33)  is  used  for 
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this  work  to  generate  the  best  FM  dorrain  possible.  Whatever 
rule  is  used,  once  the  exponential  sample  points  have  been 
corrputed  for  an  FM  preprocessor,  they  needn't  be 
recorrputed,  but  can  be  stored  for  rapid  access  during  the 
interpolation.  Soire  interpolation  rrust  be  perforrred  to 
approxirrate  the  spectral  values  of  the  new  sample  points. 
Third  or  even  forth  order  Lagrange  polynomials  have  been 
recommended  and  used  for  this  purpose  with  apparent 
success  [9,10,19].  The  advantage  of  the  Lagran.^e  technique 
is  that  its  notation  is  particularly  compact,  consisting  of 
simple  summations  and  repeated  products  that  lend 
themselves  to  digital  realization.  Unfortunately,  for  the 
data  sets  used  in  this  thesis,  the  third  order  algorithm 
was  observed  to  behave  poorly,  adding  a  ripple  in  regions 
of  rapidly  changing  slope.  That  this  might  be  the  case  was 
suggested  by  the  advice  that  Lagrange  is  very  good  near  the 
central  data  point  when  the  order  of  the  polynomial  is 
known  to  be  the  same  as  the  order  of  the  approximation, 
otherwise  it  is  best  left  alone  [20,21]  .  In  its  place,  a 
second  order  polynomial  was  used  to  interpolate  the  warped 
samples  with  results  that  were  nearly  indistinguishable 
from  the  actual  waveform,  as  seen  later  in  figure  8. 
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3.    ZEEO   POINT    CORRECTION 

Another  problem  becorres  apparent  when  the  Mellin 
transforir   of   h(f)    is   recalled, 

H a,  =   j]*^  f'Mf  -  \\c.-)  e-^*  (34 ) 

where  s=-J2f>m.  The  exponentially  sarrpled  waveforin 
described  above  is  applied  to  an  FFT  block.  As  f  approaches 
the  folding  frequency,  h(f)  tends  to  zero.  Unfortunately, 
as  f  approaches  zero,  the  value  of  h(f)  is  not  zero.  In 
fact  it  is  frequently  rather  high.  To  make  matters  even 
more  interesting,  the  left  integral  in  equation  (34^ 
clearly  shows  that  the  closer  to  zero  f  gets,  the  more 
important  h(f)  becomes  to  the  integral.  Several  solutions 
to  this  problem  are  considered  below. 

One  practical,  simple  approach  is  to  set  the  DC  (ie, 
f=0)  term  of  the  FFT  to  zero.  The  effect  is  nothing  more 
than  removing  a  DC  level  back  in  the  signal  domain,  but 
Mellin  transformation  into  the  FM  domain  leaves  the 
spectrum  dependent  upon  the  scale  factor  k  [22j .  Setting 
the  f=0  coefficient  to  zero  corresponds  to  setting  h(f)  to 
zero  for  0   f<l,  where  the  unity   upper   limit   is   chosen 
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without   loss   of   generality.  The  resulting  transform  of  a 
scaled  signal  doirain  h(f/k)  is 


^  i  (fM)f'V^=  Jc'Cjl^^^^'''jf 


(35) 


which  Is  obviously  dependent  upon  k.  The  closer  k  is  to 
unity,  the  srraller  the  effect.  By  increasing  the  f  spectral 
resolution,  the  error  can  he  reduced.  The  error  may  he 
insignificant  for  many  applications  [5-8],  "but  the 
technique  should  he  used  with  care. 

The  first  solution  has  highlighted  the  need  for  a  zero 
point  correction.  Another  common  solution  [10-11]  is 
developed  hy  breaking  the  integral  up  as  before.  Again 
using  unity  as  the  upper  limit  of  the  left  Integral  while 
remaining  general  in  application, 

i  ^(;,',o-Jt  X"^^^^  ^'"^^  (36) 

Two  assumptions  are  made  to  get  the  correction  term.  First, 
that   h(f)   remains   a  constant  h(0)  over  the  interval  of  f 


between  zero  and  one.  Second,  and  not  as  easily, accepted, 

(37) 


Equation  (37)  pretty  well   shows   why   this   assumption   is 
suspect,  but  playing  along  for  the  momrrent ,  the  question  is 
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reserved  for  a  later  detailed  look.  Accepting  the 
assurrption,  the  correction  factor  tecomes. 

Because  H(s)  is  a  corrplei  function,  Z(s)  rrust  be  applied 
(added  to  the  imaginary  part  of  the  succeeding  Fourier 
transforrr)  before  the  iragnitude  is  take  to  remove  scaling 
dependence.  This  correction  is  specifically  derived  for  use 
with  a  continuous  Fourier  transform  such  as  the  optical 
Fourier  in  imaging  systems  with  the  added  stipulation  that 
h(f)  he  nearly  constant  over  the  range  0<f<k  where  k  is  the 
largest  scale  factor  expected.  If  an  FFT  is  employed  to 
make  this  final  transform,  another  correction  should  "be 
applied  as  shown  hy  Zwicke  and  Kiss  [11]  helow.  This 
correction  factor  differs  from  the  first  due  to  an 
invariant  property  of  the  FFT.  The  FFT  of  two  unit  step 
functions  that  vary  only  in  scale  are  balanced. 


(3S) 


where  m  and  p  are  arbitrary  integers  greater  the   zero   and 

less   that  ^.      Successive  FFT  coefficients  are  sumrred,  and 

the  average  value  of  the  contributing  terms  taken  resulting 
in. 


r  ~  -t-  ^^'-  ^  "  -^  '"*'  (40) 
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This    is    then   rmjltiplled   hy  h(0)    to    arrive   at    the   FFT   ^'ellin 
correction    factor. 


•PFT 


(41) 


When  k/M  is  small  the  imaginary  term  dominates  and  the 
correction  factor  approaches  that  used  in  the  continuous 
case  (38).  Most  of  the  work  done  for  the  thesis  on  the 
exponential  algorithm  was  done  using  the  inappropriate  zero 
point  correction  (38).  Since  its  discovery  was  coincident 
with  that  of  more  powerful  methods  discussed  in  Chapter 
III,  little  data  was  taken  using  (41). 

To  hound  the  error  involed,  an  acceptable  h(t)  is 
defined,  windowed  and  transformed  using  the  Mellin 
integral.  The  window  limit  is  then  allowed  to  grow 
unbounded  and  the  resulting  expressions  are  interpreted  as 
the  error.  It's  heen  assumed  that  for  0<t<l,  h(t)=h(0). 
Since  the  error  resulting  from  the  assumption  in  equation 
(37)  arises  from  this  interval  alone,  t>l  is  ignored  for 
the  time  helng.  Warping  the  signal  as  hefore  h(x)=h(0)  for 
all  x<0  and  h(x)=0  for  all  t>0.  Evaluating  the  integral  to 
a  finite  window  width  (T) , 


/coi  /,  _  ^-^^^) 


-T 


s;^C<^r^i.)\l    -i'-^"^2. 


(42) 
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The  term  in  brackets  is  the  magnitude  of  the  contribution 
from  the  region  of  integration,  -T<i<0.  Note  that  it 
involves  a  sin{)/()  term.  The  effect  is  to  add  a  peak  of 
(T)h(0)  at  the  origin.  The  size  of  the  peak  depends 
directly  on  the  window  border  T.  Letting  T  approach 
infinity  raises  the  spike  at  u>=0»  with  lesser  peaks  at 
«^  =(2i+l )  tr/T,  where  i  is  an  integer.  Each  of  the  subpeaks 
has  a  magnitude  of  2h(0)T/(  (ei+DTT)  .  Substituting  u5  into 
the  second  relation  yields  the  envelope  (in  brackets)  and 
phase  as  T  tends  to  infinity. 


(43) 


T'-*<» 


The  error  bound  in  brackets,  does  not  depend  on  the  sample 
rate,  or  the  size  of  the  window.  Any  approximation 
approaching  zero  will  have  the  same  bound.  Although  the 
magnitude  of  the  error  is  in  a  convenient  form,  the  phase 
Is  indeterminate.  For  the  correction  to  be  applied,  the 
complex  addition  must  occur  prior  to  the  modulus  being 
taken.  This  cannot  be  done,  leaving  the  error  uncorrected, 
but  is  bounded  by  2h(0)/«O  for  continuous  and  aueriodic 
discrete  Fourier  transforms.  For  the  FFT,  equation  (43) 
does  not  bound  the  error.  The  sum  of  1/i  does  not  converge 
as  1  tends  to  infinity.  At  any  point  on  the  FFT  this  sum  is 
present  due  to  the  apparent  folding.  The  error  itself  is 
not  unbounded  because  phase  differences  in  the  sum   of   the 
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errors  at  any  point  rray  result  In  the  envelopes  adding 
destructively,  reducing  the  actual  error.  Adding  the 
envelopes  is  an  unrealistic,  worse  case  approach.  The  FFT 
correction  (41),  can  not  he  compared  to  (43)  for  these 
reasons.  The  other  error  correction,  using  the  assurr^ption 
in  equation  (37),  can  he  compared.  The  first,  setting  h(0) 
to  zero,  or  Just  plain  ignoring  the  0<t<l  interval  have  the 
error  function  hound  in  equation  (43).  Although  only 
differrlng  by  a  factor  of  two  in  magnitude,  the  constant 
phase  is  arbitrary,  and  equivalent  to  setting  T=0;  that  is, 
assuming  h(0)=0  over  2<t<l.  This  was  the  very  problem  the 
correction  was  developed  to  remedy,  tut  is  without  effect. 
Equation  (41),  the  correction  for  the  IFT  is  not  completely 
accepted  by  this  author,  albeit  no  real  empirical  evidence 
have  served  to  verify  or  dispute  the  claim.  Suspicions  are 
raised  on  two  aspects.  The  indeterminate  phase  of  the  error 
(43)  In  the  continuous  case  arises  naturally  from 
approaching  the  t=0  limit  with  the  Mellin  integral.  This 
quality  is  conspicuous  by  its  absence  in  the  FFT  error 
correction.  The  FFT  error  correction  is  computed  by  summing 
the  FFT  coefficients  in  the  complex  plane.  The  average 
position  of  the  resulting  polynomial  is  the-  error 
correction  term.  Fowever,  the  FFT  coefficients  of  a  finite 
duration  signal  are  the  values  of  the  signal,  evaluated  at 
M  evenly  spaced  points  about  the  unit   circle   [3,23].   If 
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the  sequence  is  a  constant,  the  average  value  Is  the  origin 
of  the  corrplex  plane. 

The  zero  point  corrections  for  the  Mellin  transforrr  are 
unbounded  at  t*^=0.  tiore  tirre  could  have  "been  spent 
deterrrining  the  best  applied  correction  for  the  specific 
case  at  hand,  "but  direct  methods  are  developed  in  the  next 
chapter  that  obviate  the  need  to  errploy  the  correction  at 
all. 

C.  TESTS  AND  RESULTS 

It  is  worth  admitting  at  this  point  that  the  results 
using  the  exponentially  warped  algorithm  to  achieve  a 
discrete  Mellin  transform  have  not  been  good.  More  recently- 
developed  techniques  in  Chapter  III  greatly  surpass  the 
results  reported  in  this  section.  Although  much  of  the 
theory  used  to  improve  performance  in  the  following 
chapters  could  have  been  used  here,  this  was  a  preliminary 
attempt  that  was  later  abandoned.  The  FM  processor 
described  in  the  previous  section  was  built  using  FORTRAN. 
Appendix  A  is  the  documented  program.  This  section  will 
review  the  processing  with  actual  plots  of  the  signatures 
at  different  stages,  discuss  the  required  tests,  introduce 
the  testing  approach,  and  finally  intepret  the  results.  The 
functional   block  diagram  of  the  preprocessor,  figure  5,  is 
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near  the  beginning  of  this  chapter  and  could  be  profitably 
reviewed.  Figures  6  through  9  represent  the  step  by  step 
view  of  the  signal  processing,  where  the  signal,  a  test 
shape,  is  shown  in  figure  6.  The  vaveforrr,  appears  as  an 
envelope,  and  is  drawn  with  vertical  lines  indicating  the 
sarrpled  series.  Figure  7  is  a  picture  of  the  FFT,  in  its 
sainpled  version  showing  its  symmetries.  Figure  8  is  a  very 
close  approxirration  to  the  continuous  periodic  transform 
achieved  by  filling  zeros  to  obtain  the  requisite  spectral 
resolution.  A  '*  +  "  on  the  transform  plot  indicates  an 
exponential  sample  point  interpolated  from  the  sixteen 
unique  points  in  figure  7.  The  warped  samples  are  sent 
through  the  FFT  block  once  more  with  the  result  shown  as 
figure  9.  Heavy  spectral  coloring  by  the  h(0)/i*i  correction 
factor  is  evident.  Only  the  first  half  of  the  spectrum, 
from  zero  to  the  folding  frequency,  is  valid. 

Coinplete  scale  invariance  was  never  realized  although 
its  effect  was  greatly  reduced.  All  the  testing  done  on  the 
exponential  algorithm  was  an  attempt  at  achieving  and 
verifying  shift  and  scale  invariance.  Much  effort  was  spent 
structuring  the  tests  to  avoid  the  effects  of  processing 
noise  so  the  actual  algorithmic  characteristics  could  be 
determined.  Along  the  way,  requirements  arose  to  select  a 
suitable  interpolation  polynomial,  an  optimal  zero  point 
correction  and  other  system  improvements.  In  all  cases,  the 
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Test  Fourler-Mellln   Eorrala 
Figure   9 
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test  procedure  was  to  first  select  canonic  test  shapes. 
Squares  and  triangles  were  most  frequently  used,  being 
shifted,  scaled  and  corrbined  to  determine  preprocessor 
characteristics.  Scaling  was  most  frequently  by  a  factor  of 
two  or  less.  This  corresponds  to  an  aspect  angle  change  of 
60  degrees  from  the  unsealed  case.  For  many  tests,  care  was 
specifically  taken  to  exactly  scale  a  sampled  series 
instead  of  the  waveform  envelope.  The  variation  in  input 
signal  when  this  isn't  done  is  evident  when  considering 
figure  10.  For  the  envelope  shown,  the  sampled  series 
cannot  reconstruct  the  same  signal,  and  may  result  in 
feature  space  variation.  Two  feature  qualities  were 
monitored  in  each  test  to  determine  preprocessor 
performance;  insensitlvity  to  shifting  and  scaling,  and  the 
ability  to  differentiate  between  canonic  classes.  These 
qualities  were  measured  by  visual  comparison  of  the  FM 
features,  by  computing  the  correlation  coef f iecients  and 
the  mean  squared  error  between  the  feature  sets  of 
differently  scaled  similar  test  shapes,  and  by  computing 
the  distribution  of  the  error  over  the  feature  space.  The 
latter  test  was  an  attempt  to  locate  feature  regions  of 
class  commonality,  and  regions  of  distinction  between 
canonic  classes.  This  approach  allowed  macroscopic  and 
microscopic  examination  under  changes  of  scale,  shift  and 
shape.   The   clustering   and   separation  qualities  are  data 
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dependent.  Since  no  ship  radar  video  data  was  used,  all 
observations  were  with  respect  to  canonic  classes.  No 
strong  groupings  were  detected. 

As  alluded  to  earlier,  the  results  for  the  exponential 
algorithm  were  less  than  satisfactory.  Test  shapes  were 
carefully  designed  to  irlnirrize  sarrpling  effects,  and  to 
ensure  low  side  lohes  in  the  frequency  domain.  Many  tests 
were  conducted  using  each  of  the  zero  point  correction 
methods  with  varied  shapes,  sample  rates,  and  spectral 
resolutions.  Classification  on  the  basis  of  signal  shape 
was  very  poor.  The  strongest  correlation  was  between  shapes 
of  common  duration.  Table  1  shows  a  typical  result  of 
comparing  a  rectangular  shape  and  a  raised  ramp.  The 
scaling  in  each  case  was  by  2  (60  degrees).  Equation  (36), 
Z=h(0)/ui  was  used,  but  the  others  offerred  little 
improvement.  Consistently,  the  strongest  similarity  was 
shown  between  shapes  that  had  the  same  sample  length,  vice 
shape . 
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TABLI  1 

Canonic  Shape  Fourier  -  Mellin 
Feature  Corrnarisons 


a.  Peak  Correlation  Values 


RFCT 

RECT/2 

RAMP 

RAMP/2 

RECT 

1.00 

0.77 

0.98 

0.76 

RECT/2 

- 

1.00 

0.87 

1.00 

RAMP 

- 

- 

1.00 

0.86 

RAMP/2 

— 

— 

— 

1.00 

b.  Squared  difference  "between  features. 


RECT 

RSCT/2 

RAMP 

RAMP/2 

RECT 

.000 

.032 

.019 

.032 

RECT/2 

- 

.000 

.032 

.019 

RAMP 

- 

- 

.000 

.032 

RAMP/2 

— 

_ 

— 

.000 
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Ill .   DIRECT  MELLIN  TRANSFORMS 

The  last  chapter  developed  a  rrethod  of  obtaining  the 
Mellin  transform  hy  exponentially  warping  the  signal  prior 
to  using  an  FFT  hlock.  This  technique  is  referred  to  as  the 
fast  Mellin  transform  (FMT).  Although  the  promise  of  scale 
invariant  features  is  attractive,  some  of  the  prohlems 
encountered  that  make  the  FMT  unattractive  are  reviewed 
here.  The  required  sample  rate  varies  with  respect  to  the 
data,  making  general  applications  difficult.  The  tendency 
is  to  use  more  samples  than  required,  which  quickly  hecorres 
costly  in  an  exponential  sampling  scheme.  The  need  to 
exponentially  warp  (to  interpolate)  a  set  of  new  samples, 
is  expensive  in  time  required.  For  true  scale  iavariance,  a 
correction  factor  is  required  hut  because  of  the  integral's 
unbounded  nature  at  zero  this  correction  factor  is 
indeterminant.  Several  correction  methods  have  been 
employed,  but  they  depend  on  unspecified  data 
characteristics.  These  effects  combine  to  irake  the  actual 
performance  of  the  algorithm  poor.  Although  scaling  effects 
are  mitigated,  they  remain  an  artifact,  which  is  disturbing 
to  classification  attempts.  This  chapter  outlines  the 
effort  to  remove  these  limitations.  Some  useful  Mellin 
properties  are  developed,  and   then   applied   to   establish 
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several   Direct   Mellin  Transforms  (IMTs)  which  were  built, 
and  their  performance  compared. 

A.  SO^'E  USEFUL  PROPERTIES 

Sorre  general  observations  are  made  here  about  the  Mellin 
transforms,  and  are  followed  by  some  specific  relationships 
which  were  derived  and  applied.  A  property  of  the  Mellin 
s-domain  is  that  it  is  unaffected  by  scaling  changes  in  the 
original  x-domain.  Figure  11a  is  a  test  shape  in  the 
x-domain.  Two  features  are  identified  according  to  their 
amplitudes  A  and  B,  at  x=a  and  x=b  respectively.  The  ratio 
of  a/b  equals  c.  To  be  simply  scaled  by  k,  h(x/k)  must 
remain  the  same  in  all  aspects  except  that  the  distance 
between  features  has  been  changed  according  to  the  scale 
factor  k.  Figure  lib  shows  the  scaled  domain.  Features  A 
and  B  are  again  identified  at  their  scaled  positions  ka  and 
kb.  The  signal  property  maintained  by  scaling  is  relative 
positional  integrity,  that  is,  ka/kb  equals  c  as  before. 
The  positional  integrity  of  the  features,  the  ratio  of 
their  distance  from  the  origin,  to  that  of  another's  is 
unchanged.  Restated,  an  operation  0[h(x)l  in  the  x-domain, 
will  leave  the  s-domain  modulus  invariant  to  simple  scaling 
by  k  if 

0[iC^/A)]'=    f('x/A)  (44) 
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Note  that  the  entire  dorrain  is  scaled,  so  no  unsealed  shift 
in  the  dorraln  can  he  permitted  as  already  discussed.  In  the 
Mellin  s-dorrain,  any  simple  scaling  in  the  x-domaln  results 
in  a  phase  distortion  (the  modulus  is  invariant  to  k). 
Manipulations  in  the  s-domain  will  leave  that  domain 
invariant  to  k,  as  long  as  the  modulus  is  modified  hy  a 
multiplicative  factor  of  constant  phase.  That  is,  if  !G(s)| 
is  an  arhitrary  function  of  s  (except  that  it  does  not 
depend  on  k)  and  [M  [h(x/k)]  j  =  !H( s) i  is  also  invariant  to  k, 
then  their  product  is  invariant  to  k  and  the  x-domain 
remains  simply  scaled.  For  instance,  in  the  x-domain  the 
operator   0[h(x)] =x(h(x) ) ,    does   not  meet    condition    (44). 

^A(,yxM)      7^     /"C-^Z-^J  (45) 

So  the  l^ellln  transform's  modulus  of  xh(x/k)  cannot  te 
invariant  to  k. 

In  Chapter  II,  a  error  was  apparent  hecause  h(x)  was  not 
equal  to  zero  for  x=0.  If  h(x)  could  he  modified  with  an 
operation   that   met   the   condition   of   equation  (44)  and 
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always  produced  a  series  that  was  zero   at   x=0   a   general 
approach  can  be  developed.  Consider  two  operators. 


0[1C^/A)]    =  ^C-^Ca/^})=  /C^/^) 


(46) 


J 


Ol-AC'K/A)]^  ^  j:;(.k(:^/A))^  fc^/A) 


(47) 


Equation  (46)  will  produce  an  acceptable  f(i)  as  df/dx=0  at 
x=0.  Equation  (47)  will  always  produce  the  required 
condition.  To  see  the  frequency  dorrain  equivalents,  we  rrust 
assuire  that  the  Mellin  integral  exists.  Integrating  by 
parts , 


\   \  Jl(^/A)  ^^'^d 


-    (//of 


4^ 


I [  ni'J.<:^/A)] '-  s  j"^ U/4 ) ^ ^"^ 


==j'S  ^^M^/A)  ^'-U^  \^\'  ^^'^1 


(48) 


The  result  in  the  frequency  domain  is  consistant  with  the 
conditions  stated  above.  The  limit  of  h(x)  as  x  goes  to 
zero   or  to  infinity  must  be  zero  if  equation  (48)  is  to  be 
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true.  Similarly,  under  the  sarre  conditions,  it  can  be  shown 
that 


l'>nli;c^ic..A))]\=  /^-^-^  ^'^>\ 


(49) 


Equations  (48)  and  (49)  clearly  do  not  apply  where  h(0)  is 
not  equal  to  zero.  However,  "by  using  the  x(  d(  h(x)  )/dx) 
irodlfying  operation  of  (47),  a  series  can  always  "be  zero  at 
1=0.  The  function  h(x)  is  further  constrained  "by  the  fact 
that  it  must  fall  off  faster  than  1/x.  This  assumption  must 
be  valid  and  the  modifying  operator  must  he  applied  in  the 
x-domain  for  the  Mellin  integral  to  exist  in  general.  A 
Mellin  transform  of  a  function,  after  having  a  modifier 
applied,  will  he  called  a  modified  Mellin  transform. 


dnc        ^   ^^  (50) 

The  integral  is  close  to  a  form  which  is  realizable,  except 
for  the  upper  limit.  For  a  finite  sampled  series,  h(n)  will 
he  assumed  zero  outside  of  the  interval  0^n<N.  This 
truncation  effects  the  transform  of  an  otherwise  infinite 
series.  In  this  application  the  Mellin  is  applied  to  an  FIT 
frequency  spectrum.  The  truncation  in  the  frequency  domain 
is  due  to  hand  limiting  the  signal  prior  to  sampling. 
Scaling  in  the  time  domain  will  not  result  in  simple 
scaling,  hut  will  add  a  dependence  on  the  scale  factor  k 
that  can  not  he  removed  by  the  transform.   An  approach  by 
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Prost  and  Goutte  is  used  to  predict  the  size  of  the  error 
[24,25].  First  a  suitable  function  will  be  selected  and  a 
relative  error  of  truncation  (RET)  deternrined  and  applied 
to  two  scalings.  The  relative  difference  of  the  feature 
space  is  found  and  Identified  as  the  error.  Remerrhering 
that  this  is  applied  to  a  frequency  spectrum,  dh(f)/df  is 
approximated  by  a  function  of  the  form. 

J.  1(4)    ^    r  --r 
^7?  ^  ^  (51) 

The  modified  Mellin  transform  of  (51)  over  a  finite 
frequency  range  would  be  approximately. 


H 


4. 


(52) 


The  lower  limit  has  been  set  in  a  manner  to  be  consistant 
with  Plancherel's  theorem.  A  convenient  worse  case 
assumption  is  that  the  lower  limit  is  essentially  zero,  but 
this  depends,  in  general,  on  7  and  the  data  itself. 
Ha*(s,f)  converges  toward  Ha ( s )  ,  equation  (50),  in  the  mean 
square  as  F  tends  to  infinity.  The  mean  square  error  will 
not  provide  an  expression  for  the  error  that  can  be  used  to 
correct  for  it  but  is  useful  to  measure  the  effect  of  the 
truncation  in  more  general  terms.  A  relative  error  of 
truncation  (RET)  can  be  computed  if  the  error  is  assumed  to 
be  distributed  evenly  over  the  range  from  which  Ha*(s,f)  is 
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computed.   By   employing  Plancherel 's   theorem,   the   mean 
squared  values  are, 


^o  -J  a 


(53) 


=/^  rr=-=/^ 


(54) 


The  RET  is  defined  and  solved  as. 


K^r  ^    I  ^l^'  ]  ^=  g''^(zf\F  ^  l)  ^  (55) 


But  the  limit  F  depends  not  only  on  the  pass  band  F,  hut  on 
the  scaling  in  the  f-domain.  A  relative  error  (R2)  between 
a  truncated  spectrum  and  a  scaled  and  truncated  signal 
would  he  more  complex,  hut  worth  the  effort. 


where  0<k<l .  Two  observations  should  be  made  here.  First, 
the  relative  error  of  truncation  (55)  and  the  relative 
error  or  difference  between  two  truncated  spectrums  scaled 
differently  (56)  both  depend  on  F.  F  depends  on  the  cut  off 
frequency  of  the  low  pass  filter  prior  to  any  sampling.  It 
is   often   chosen   to   minimize  aliasing  depending  upon  the 
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limitations  of  the  sampling  circuit.  Second,  RI  depends  on 
k  as  well.  The  importance  of  scaling  differences  to  this 
data  type  can  be  readily  seen.  If  equations  (55)  and  (56) 
are  valid,  and  if  the  range  of  k  can  he  hounded  (a  design 
specification)  F  can  he  chosen  to  realize  a  stated  RE.  Or 
if  F  is  fixed  and  k  hounded,  the  RE  may  he  determined  for 
evaluation.  If  the  data  type  is  not  appropriate,  a  more 
representative  function  may  be  determined  and  used  in  place 
of  equation  (51)  to  attain  better  expressions  for  RET  and 
RE. 

3.  ALGORITHM  DEVELOPMENT 

Part  A  above  provided  some  background  for  making  the 
Eirect  Mellin  Transforms  (EMTs)  in  this  section.  In  this 
presentation  the  Justification  is  given  with  the 
application  first.  However,  chronologically  the  neat 
package  presented  above  was  preceded  by  extensive 
evaluation  of  empirical  results,  not  vise  versa.  Appendix  3 
documents  the  FORTRAN  implementation  of  all  the  algorithms 
developed  in  parts  1  and  2  below.  Figure  12  is  a  functional 
block  diagram  of  a  generalized  FM  preprocessor  using  a  DMT 
to  show  the  simplicity  with  which  it  can  be  applied,  as 
opposed  to  the  FMT  covered  in  the  second  chapter,  figure  5. 
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1 .   First  Difference  Approximations 

Although  developed  through  a  different  rationale, 
this  first  algorithrr  was  developed  ty  Zwicke  and  Kiss  [ll]  . 
Starting  with  a  sampled  series  hi,  the  series  is  operated 
on  by  the  rrodifier  defined  in  equation  (51),  using  the 
first  backward  difference  to  approximate  the  derivative 
with  respect  to  x.  Unit  step  size  is  assumed. 

^    jA->  =    ^  (i^-A^J   =  ^  VJ^  (57) 

Taking  the  trapezoidal  rule  to  evaluate  the  modified  Mellin 
integral  (52)  while  recalling  that  h(0)=0,  and  h(N)  is 
assumed  zero, 


H.,  <^--i=  E  vA^^^ 


where  s=~j2'f^m/M.  The  complex  coefficients  are 


(58) 


/n 


'=i  cos  Z'ff'  C'»^  /m  )  ^  ^  "d  -^"^  Z<f  C'^  /m  )A^  C^) 


(59) 


They  can  "be  calculated  off  line  and  stored  to  produce  just 
the  desired  characteristics.  The  factor  (2''^m/M)  could  be 
any  number  that  produces  an  interesting  feature.  Undesired 
features  need  not  be  computed  (what  in  general  was  a  N  by  f^ 
matrix,  where  M  is  the  number  of  Mellin  transform 
coefficients  and  N  the  number  of  x  sample  points).  If  a 
relatively  small  number  of  features   is   required,   perhaps 
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the  processing  will  be  manageable.  Notice  that  no  zero 
point  correction  is  required.  Cnly  data  changes  contribute 
to  the  transform.  These  observations  are  valid  for  any  of 
the  rrodified  Direct  Mellin  Transforrrs  developed  in  this 
section. 

By  using  a  central  difference  instead  of  the 
backward  difference,  a  siirlliar  result  is  obtained. 

H..C^)  =.  ^  U^^^  -^^..  )  -^  V2  (60) 

Other  nuirerical  integrations  tnay  be  used  with  irrproved 
results,  and  other  rrethods  can  be  used  to  increase  the 
order  of  the  approiirration. 

To  test  the  algorithms,  a  ranp  and  inverse  ramp  were 
used  in  a  scaled  and  unsealed  rr^ode.  The  rarrps  and  their 
scalings  are  shown  in  figure  13.  Figure  14  is  the  analytic 
results  of  both  waveforms  plotted  with  the  transform  found 
using  equation  (60).  This  is  a  dramatic  improvement  over 
anything  used  with  the  methods  discussed  in  the  previous 
chapter.  The  noise  of  the  signal  appears  to  be  diverging  as 
frequency  grows,  but  over  the  range  plotted,  the  appearance 
is  that  of  an  algorithm  trying  to  do  well. 
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2 .      Second  lifference  Approximations 

A  second  difference  algorithn^  can  be  achieved  by- 
proceeding  as  before.  The  pure  scaling  operator  used  to 
prepare  the  signal  is, 


OC^C^jJr  ^'^^    =  ^^(^^..-Z^^^^KJ 


(61) 


And  the  new  algorithm  is. 


H 


*.z 


S-^l 


(62) 


Other  second  order  operators  have  been  used,  and 
partitioned  into  forward  and  backward  difference  variations 
to  (61),  but  this  appears  to  be  a  basic  and  useful  form. 
The  color  l/(s-»-l)  is  present  to  approximate  the  modified 
Mellin  of  equation  (60).  The  term  (s+1)  is  valid  assuming 
that  dh(x)/dx  is  exclusively  upper  bounded  by  ln(x)/x  as  it 
approaches  zero  or  infinity.  For  comparison,  another  second 
difference  algorithm  was  developed  based  on  the  modifier 
(x(d/dx)(x(d/di)h(x). 

This  is  roughly  the  sum  of  the  methods  defined  in  equations 
(60)  and  (62)  above.  The  assumption  for  deriving  the  color 
1/s  is  even  less  restrictive  than  before,  but  the  algorithm 
performs  poorly  compared  to  (62).  These  transforms  depend 
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on  second  difference  characteristics,  but  are  modified  by  a 
1/s  terrp  which  has  a  stabilizing  effect.  This  is  equivalent 
to  a  division  by  x  and  integration  in  the  x-domain .  The 
order  of  the  approximation  has  been  increased  by  the 
rrodifier.  Results  using  equation  (60)  and  (62)  should  be 
alike.  Figure  15  is  the  results  of  using  (62)  cofrpared  to 
the  closed  form  solution  to  the  figure  13  test  shape 
transforms.  An  improved  performance  over  (60)  is  seen  over 
some  of  the  range,  but  a  drop  as  frequency  increases 
degrades  the  accuracy  in  figure  15b. 

3.   Higher  lifference  Approximations 

Higher  difference  approximations  can  be  developed. 
For  instance,  one  algorithm  depending  upon  the  third 
difference  is, 

(64) 


H.3^^£^'^--"^ 


The  performance  of  higher  order  algorithms  becomes 
Increasingly  suspect  because  of  the  extreme  weighting  they 
apply  to  different  parts  of  the  series.  Different 
algorithms  exist,  but  this  weighting  is  always  a  factor.  A 
smoother  transform  is  achieved,  but  a  large  error  is  likely 
to  develop  due  to  the  algorithm's  dependence  on  higher 
order  derivatives  and  the  nature  of  sampled  data.   However, 
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300 


these   corrrrents    are   speculative   since   they   were   not 
coafirrred  experimentally. 

4 .   Higher  Order  Integrations 

Higher  order  integration  rules  should  "be  atle  to  te 
used  with  a  corresponding  improvement  in  performance.  One 
using  the  first  difference  with  Simpson's  rule  was 
implemented  with  di ssappointing  results.  Figure  16  is  the 
result  of  such  an  implementation.  The  droop  for  the  ramp 
input  is  apparent  even  though  not  present  in  the 
trapezoidal  rule  used  in  subsection  1  above.  The  higher 
frequency  error  is  also  more  prevalent  than  before.  A 
program  error  is  of  course  suspected,  but  was  never  found. 
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IV.  CLASSIFICATION  PREPROCESSING 

The  problem  of  determining  important  signal 
characterisics  can  be  approached  in  at  least  tvo  different 
ways.  First,  by  trying  to  learn  what  is  important  in  human 
recognition,  and  then  trying  to  adapt  a  machine  to  emulate 
that  behavior.  Or  second,  by  using  successive  transforms  to 
remove  information  known  to  be  superfluous  to 
classification,  while  keeping  enough  information  to 
reliably  assign  an  object  to  a  class.  Addressing  the 
former,  even  though  it  is  difficult  to  determine  specific 
details,  some  key  aspects  of  human  visual  recognition  are 
discernable.  Chief  among  these  is  that  the  intensity  level 
of  a  scene,  or  object,  does  not  appear  to  be  as  important 
as  the  relative  position  of  the  edges,  i.e.  the  shape 
separating  different  intensities  and  frequencies  [25-28]. 
Examples  in  scene  analysis  show  clearly  that  the  edges  or 
shapes  are  far  more  critical  to  human  recognition  than  the 
relative  power  differences  themselves.  The  invariant  shapes 
or  angles  in  scene  analysis  find  their  analog  in  ratios 
between  similiar  points  in  differently  scaled  time  series. 
Examination  of  many  range  only  radar  video  ship  signatures 
has  provided  a  basis  for  noting  that  the  relative  position 
of  time  domain  features  remain  constant  over  changes  in 
aspect   angle,  while  the  relative  intensity  of  the  features 
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vary  greatly  as  shown  in  figure  17.  As  the  aspect  changes, 
if  this  set  of  ratios  remains  constant  within  a  ship  class, 
then  the  set  is  identified  as  inforrration  worth  preserving 
for  the  classifier.  Conversely,  since  relative  intensitiy 
is  not  a  stable  measure  over  aspect  angle,  or  from  among 
different  ships  of  one  class,  that  information  should  te 
intentionally  removed  to  provide  tighter  natural 
clustering,  with  the  minimum  number  of  features. 

A.  INFORMATION  REQUIREE  TO  CLASSIFY 

There  are  two  preconditions  that  must  both  exist  for  a 
set  of  possible  input  signatures  to  separate  into  distinct 
classes.  First,  the  features  of  a  particular  class  must 
have  some  common  characteristic  about  them,  and  second, 
this  characteristic  must  in  some  way  be  unique  with  respect 
to  other  classes.  The  assumption  in  existing  radar 
signature  classification  projects  is  that  there  is  enough 
Information  in  the  signatures  to  permit  this 
classification.  Short  of  actually  trying  to  classify  with  a 
set  of  realistic  signatures,  the  analytical  determination 
that  sufficient  information  is  present  in  a  set  of  all 
possible  ship  signatures  is  difficult  to  approach. 

To  establish  how  well  the  autonomous  classifier  is 
performing  some  measure  of  the  classif iablli ty  of   the   set 
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received  signals  is  desired,  failing  this,  some  discussion 
of  the  inforrration  capacity  of  the  preprocessor  should  at 
least  be  considered.  The  preprocessor  produces  the  FFT 
rragnitude  of  a  sarrpled  signal  as  the  output  of  the  first 
stage.  Most  of  the  unique  positional  relationships  of  the 
signal  upon  which  hurran  recognition  apparently  depend  has 
been  destroyed.  Next,  the  Mellin  transfrom  stage 
effectively  distorts  the  signal  and  uses  the  rragnitude  of  a 
second  Fourier  transform  as  the  output  features.  Signals 
reconstructed  on  the  basis  of  FFT  phase  information  alone 
usually  provide  sufficient  similarity  to  be  associated  with 
the  original  signal,  whereas  reconstruction  on  the  basis  of 
magnitude  does  not  retain  any  significant  detail  except 
when  the  signal  is  symmetric  [29] .  The  data  is  also  masked. 
For  most  applications,  the  data  is  frequency  shifted, 
filtered  and  sampled  as  a  baseband  signal.  This  windowing 
masks  the  magnitude  characteristic  and  is  specifically 
designed  into  the  processing.  The  resulting  FM  features  are 
insensitive  to  positional  and  scaling  relationships  that 
are  necessary  in  human  visual  recognition. 

Analytic  support  is  also  available  to  quantify  the 
importance  of  relative  position  of  events  [29].  By 
considering  rms  error  (due  to  spectral  phase  and  amplitude 
quantization  for  random  signals)  it  has  been  concluded  that 
approximately   two   more   bits   are    required    for    the 
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quantization  of  phase  inf ortration  than  arrplitude  for  the 
SAive  rrrs  error.  A  separate  analysis  applied  distortion  rate 
theory  to  real-part,  iiraginary-  part,  and  magnitude-phase 
encoding  of  the  DFT  of  random  sequences.  The  result  was 
that  phase  required  1.4  hits  iTore  storage  than  magnitude 
for  a  similar  error  [30].  A  third  approach  concluded  that 
the  Fourier  phase  includes  1.8  hits  more  information  than 
the  magnitude  [31] .  This  was  based  on  analysis  of  image 
reconstruction  from  kinoforms  (phase-only  holograms).  The 
fact  that  phase-only  reconstruction  preserves  much  of  the 
correlation  between  signals  would  suggest  that  the  location 
of  events  tends  to  he  preserved.  Further,  it  seems  that 
this  information  is  lost  by  taking  only  the  magnitude  of 
the  Fourier  transform.  Another  interesting,  albeit 
informal,  view  is  apparent  as  one  considers  the  phase-only 
signal  as  a  spectral  whitening  process. 


Ji^ct)]^  F(^)   =  \f(0\  e"? 


-.•^^ 


(65) 


For   reconstruction   by   phase  alone,    where    the     magnitude      is 
set    to   one; 


^Lf^l 
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iFC-f)l 


Jf/^o] 


(66) 


Since  the  received  radar  signature  will  have  an  abundance 
of  low  frequency  spectral  lines  and  smaller  high  frequency 
components,   the  low  frequency  information  is  not  weighted 
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as  heavily  as  the  higher  frequency  inf orir.atioa  ia  the  phase 
reconstructed  signal.  This  seems  like  it  would  accentuate 
sharp  changes  in  the  reconstructed  signature  without 
rerrovlng  relative  location  inforrration.  The  result  is  the 
surrrration  of  the  different  frequency  components,  all  with 
zero  phase  (i.e.,  no  positional  or  amplitude  distrihution 
information  can  remain). 

Considering  the  Mellin  transform  next,  in  a  continuous 
case,  the  exponential  warp  does  not  lose  any  information. 
To  zero  the  first  data  sample  as  required  by  the  Chapter 
III  modifications,  will  surely  destroy  information,  but 
this  may  be  confined  to  the  EC  term  alone.  The  information 
lost  during  the  final  transform  and  magnitude  is  difficult 
to  assess.  Increased  masking  occurs  due  to  the  spectral 
truncation,  so  actual  information  loss  may  not  be  as  great, 
but  masking  distortion  may  be  greater  than  before.  So 
approximately  two  bits  of  information  are  lost.  Only  a 
quarter  of  what  was,  remains.  Information  is  also  lost  when 
the  transforms  are  normalized  in  the  processing  so  that  any 
power  calculation  is  also  meaningless.  Some  interesting 
questions  arise.  After  removing  positional  relationships, 
scaling,  and  power,  what  signal  qualities  remain  and  are 
they  useful  in  classification?  Although  it  is  true  that 
this  insensi tivi ty  may  add   a   certain   robustness   to   the 


77 


systeir,    the   arbitrary   loss   of   valid   classification 
inf orrration  should  he  rriniFrized. 

By  exarrining  the  quality  that  rrust  be  ignored,  and  by 
comparing  its  removal  to  what  is  actually  removed  by  the 
processing,  an  interesting  result  will  develop.  The  effect 
of  a  time  domain  shift  on  the  frequency  domain  is  an 
additive  phase  term,  linearly  related  to  the  frequncy  of 
the  coefficient,  as  in  equation  (2).  Most  of  the  structure 
of  the  signal  is  held  in  the  phase  relationships  with 
respect  to  the  fundamental  and  higher  frequency  terms.  So 
shift  can  be  defined  as  the  phase  of  the  fundamental 
complex  coefficient.  3y  setting  the  phase  to  zero,  and 
adjusting  the  other  coefficients  according  to  their 
component  frequency,  the  structure  of  the  signal  is  not 
lost  but  reconstructed  about  the  fundamental  as  before.  If 
there  are  N/2  spectral  phase  angles,  only  the  fundamental 
needs  to  be  zeroed  to  remove  the  shift.  If  the  information 
is  contained  uniformly  in  the  spectral  phase,  then  only  2/N 
of  this  structural  information  needs  to  be  removed.  When 
the  magnitude  is  taken  to  produce  shift  invarient  features, 
all  of  the  phase  relationships  are  destroyed.  (N-2)/N  of 
the  information  once  held  in  the  phase  was  removed 
needlessly.  The  amount  of  information  lost  removing  the 
shift  can  be  made  arbitrarily  small.  As  N  grows  unbounded, 
the  amount  of  information  that  needs  to  be  removed  tends  to 
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zero.  This  result  seems  to  be  supported  by  experience.  With 
enough  samples  of  a  shape,  its  position  with  respect  to  the 
observation  field  is  irrmaterial.  In  theory  the  principle  is 
sound,  but  some  practical  limitations  may  degrade  predicted 
performance.  Recalling  that  an  exponential  warp  translates 
scaling  to  shifting  generalizes  the  result  a  bit  further. 


J,(t/A)  =  I  ce^-^^; 


(67) 


The  same  priciple  that  permits  simple  shift  removal  is  also 
valid  for  the  removal  of  scaling  dependence  as  well.  Using 
the  Mellin  transform,  scaling  dependence  may  be  removed  by 
zeroing  the  fundamental  and  adjusting  all  the  other 
coefficients  as  described  above.  For  the  FM  preprocessor, 
the  information  lost  removing  the  scaling  and  shifting 
dependence  may  be  made  arbitrarily  small  by  increasing  the 
number  of  spectral  samples  used.  Since  the  number  of 
spectral  samples  can  be  increased  by  filling  zeros  onto  the 
finite  signal  in  the  original  domain,  this  process  does  not 
effect  the  data  sample  rate.  This  approach  was  not  verified 
experimentally,  but  represents  a  potentially  powerful  tool 
to  analyse  and  improve  the  extracted  feature  space. 
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E.  RANGE  CNIY  RADAR 

This  section  addresses  ship  classification  on  the  basis 
of  inforiration  gathered  frorr  a  range  only  radar  video  ship 
signature.  An  exarrple  of  such  a  signature  has  already  heen 
considered  as  figure  16.  Classification  by  range  only  radar 
signatures  is  subject  to  the  same  distortions  discussed 
above.  Typically,  the  radar  return  is  detected  and  isolated 
ia  a  range  gate  that  is  sampled  and  digitized.  The 
rectangular  sampling  window  can  be  considered  the  range 
gate  itself.  The  range  gate  is  designed  to  ensure  that  the 
included  range  is  greater  than  the  maximum  ship  length  so 
that  the  time/range  windowing  has  no  effect  on  the 
frequency  spectrum  of  the  signature,  other  than  increased 
spectral  resolution.  The  placement  of  the  ship  signature  in 
the  window  is  not  set,  neither  from  encounter  to  encounter, 
nor  from  pulse  to  pulse  (Jitter).  The  total  effect  Joins 
together  to  produce  the  framing  distortions.  Signature 
scaling  results  from  viewing  the  ship  from  different  aspect 
angles.  The  sampling  rate  must  be  done  at  more  that  twice 
the  inverse  of  the  resolution  of  the  receiver.  Quantization 
levels  are  chosen  in  a  manner  to  reduce  that  predictable 
random  noise  to  an  acceptable  level.  The  pulse  to  pulse 
Jitter  is  an  ever  present  characteristic  of  the  radar 
problem,  but  at  a  normal  resolution  (greater  than  2t  feet) 
integrating  the  return   in   the   time  domain   removes   the 
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jitter  effect,  reduces  scintillation,  and  improves  the 
resolution  of  the  signature.  Although  predetection 
(coherent)  integration  is  more  efficient,  post  detection 
(noncoherent)  integration  is  rrore  corrrron  because  of  the 
convenience  of  not  having  to  preserve  the  radar  frequency 
in?)  phase.  For  post  detection  integration  of  n  pulses,  the 
signal  to  noise  ratio  would  be  something  less  than  n  times 
the  signal  to  noise  ratio  for  one  pulse  [32].  More 
important  to  the  recognition  problem  itself,  for  a  stable 
system  by  the  law  of  large  numbers  [33],  fluctuation  of  the 
average  value  of  the  return  will  be  overcome.  That  is,  with 
the  integration  of  n  pulses  the  resolution  (R)  will  become 
finer  as 


K  =  -j  S^/D 


^  (68) 


with  a  probability  of  1-E ,  where  S  squared  is  the  variance 
of  the  signal  from  pulse  to  pulse.  For  very  high 
resolution,  the  cost  of  making  the  signature  stable  with 
respect  to  the  integrator  becomes  prohibitive.  It  has 
become  convenient  to  integrate  the  spectrum  of  the 
signature  because  the  Jitter  effects  can  be  completely 
removed.  Another  limiting  factor  in  integrating  over  a 
period  of  time  is  that  the  position  and  aspect  cf  the 
target  are  dynamic.  They  change  with  time.  A  conceptually 
attractive  solution  is  a  recursive  filter  which  weights  the 
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Integrated  pulses  such  that  the  older  they  becorre,  the  less 
weight  is  accorded  to  them. 

If  the  course  and  speed  of  the  target  ship  are  known 
from  measurement  of  the  target  trade,  it  is  possible  to 
infer  the  aspect  of  the  ship.  The  range  profile  can  give 
some  estimation  of  size.  Unfortunately,  the  three 
dimensional  change  in  aspect  angle,  commonly  suffered  by  a 
ship  presents  more  than  Just  a  video  signature  scaling 
change.  The  radar  cross  section  of  even  individual 
structural  components  of  the  radar  target  changes  with 
respect  to  the  aspect  angle.  The  composite  effect  is  that 
ship  signatures  vary  greatly  with  aspect  angle.  The  radar 
is  an  electromagnetic  sensor,  reacting  to  energy  reflected 
from  the  target.  These  reflections  are  a  result  of 
scatterers  that  are  related  in  dimension  to  the  wavelength 
of  the  illuminating  energy.  Because  of  the  great  difference 
in  wavelength  between  light  and  microwaves,  what  can  be 
"seen"  by  radar  may  be  quite  different  than  that  seen  by  an 
eye.  Also,  when  measuring  size  or  any  distances  with  a 
radar  of  high  resolution  (less  than  50  feet),  an  error  can 
be  incurred  since  the  extremities  of  the  target  are  not 
always  good  scatterers.  Echoes  from  the  forward  or  stern 
portions  of  the  target  might  be  observed  in  the  noise, 
especially  for  a  relatively  low  power  radar  [32].  After 
reviewing   hundreds   of   signatures,   it  appears  that  major 
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features,  such  as  overall  ship  length  and  dorrinant  mast 
structure  are  frequently  discernable,  but  vary  in  relative 
arpplitude.  Resonance,  shadowing  (one  reflector  hiding 
another),  multipath  returns,  the  irapping  of  three 
dimensional  aspect  changes  onto  a  one  dimensional  time 
series,  and  the  amplitude  and  phase  of  component  returns 
summing  constructively  or  destructively  to  cause  a 
scintillation  of  the  composite  target.  Some  of  the 
variations  caused  by  these  conditions  can  be  lessened  by 
integration  but  major  effects  remain  causing  the  signature 
to  vary  in  shape  and  content  with  the  aspect  angle.  For 
this  reason,  a  class  feature  volume  cannot  be  reduced  to  a 
single  point,  but  will  remain  a  hypervolume  in  the  feature 
space  even  in  the  ideal  case.  Any  selection  of  features 
should  try  to  minimize  this  volume.  Features  should  be 
selected  that  are  relatively  insensitive  to  known 
superfluous  effects. 

C.  CLASS  DISCRIMINATION 

In  the  last  chapter,  the  major  concern  was  removing  two 
sources  of  variance  with  no  classifying  value.  Algorithms 
were  developed  and  canonic  shapes  generated  to  verify  the 
algorithms  and  demonstrate  the  invariance  to  scaling.  In 
the  same  manner,  this  chapter  has  reviewed  information  loss 
and  process  masking.  The  question   of   whether   sufficient 
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inf orrration  is  present  to  classify  was  raised.  To  answer 
this  question  some  sirrple  classes  of  canonic  figures  are 
defined,  and  put  through  the  entire  FM  preprocessor 
documented  in  Appendix  E.  The  LMT  algorithrr  found  to  offer 
the  best  scale  factor  rejection  in  Chapter  III  was  used  to 
generate  the  final  features.  The  algorithm  chosen  is  based 
on  a  second  difference  modification  and  is  defined  in 
equation  (62), 


H^.  C 


--T^Z^^-^^ 


/>\ 


J+l 


(62) 


'^-1 


After  canonic  tests  were  made,  preprocessor  performance  on 
several  ship  signatures  was  recorded.  Although  this  was 
premature  in  the  logical  test  sequence,  the  results  are  of 
some  interest. 

1 .   Test  Shapes  and  Results 

Four  test  shapes  were  used.  Figure  18  shows  the  test 
shapes.  All  were  scaled  and  shifted  originally  to  test  for 
algorithm  verification  and  demonstrate  scale  invariance.  In 
this  series  of  tests  they  were  left  fixed  and  used  in 
different  combinations  to  try  to  detect  shape  presence  in 
the  FM  feature  space.  The  object  is  to  differentiate 
between  different  canonic  classes.  Figure  19  shows  a 
comparison  of  the  shapes,  rectangle  and  triangle.  A  test 
combining   the   rectangle   and   triangle  was  the  subject  of 
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a.    R£Cr 


t.      TRI 


c.  SHAP£  1 


e.    SNA  PC  2 


Test  Shapes 
Figure  18 
85 


u 
111 

QC 


1.0 


0.8 


0.6 


0.4 


0.2 


0.0         I   I   I.J.  I   I   I 
0.  50. 


Mil 


Mm  I     ,l.,l      I 


too. 


150. 


200 


250. 


300 


S-MELLIN    FREQUENCY 


X 

K 


1.0 


0.8 


0.6 


0.4 


0.2 


0.0 


rMV^^TMrri  M.M 


0,  50.  100.         150.        200.        250. 

S-MELLIN    FREQUENCY 
RECT   and  TRI    IM   Features 
Figure   19 
86 


300 


figure  20.  Finally  in  figure  21  a  rectangle  with  two 
triangles  is  shown  in  the  IM  feature  space.  It  is  clear 
frorr  the  plots  that  irost  of  the  scale  variance  has  been 
removed.  Just  as  important,  some  quality  does  remain  that 
differentiates  hetween  the  canonic  shapes.  A  square  in 
general  can  he  differentiated  from  a  triangle.  A  "ship" 
with  a  single  mass  of  superstructure  can  he  separated  from 
one  with  two  such  masses. 

Although  it's  clear  from  the  plots  that  there  is  a 
unique  quality  left  in  the  feature  space  to  allow  the  time 
domain  shapes  to  he  classified,  some  quantified  measure  of 
system  performance  is  required.  The  magnitude  of  the 
feeature  vectors  have  all  been  normalized  with  respect  to 
the  first  coefficient,  so  in  that  region  little 
discrimination  can  be  expected.  For  higher  N'ellin 
frequencies,  noise  dominates.  A  region  of  coefficients, 
11-100  was  chosen  to  classify  the  shapes  by  correlation. 
The  results  are  included  as  Table  2.  The  improvement  in 
performance  over  that  shown  by  Table  1  in  Chapter  II  is 
dramatic.  The  methods  in  that  earlier  test  resulted  in  the 
observed  length  of  an  object  being  the  distinguishing 
criterion  for  classification.  Using  methods  supported  in 
Chapter  III,  variations  due  to  scaling  and  shifting  of  the 
original  domain  have  been  removed.  The  features  now  reflect 
the  shape  of  the  object  in   the   time   domain.   An   unusual 


87 


1.0 


0.8 


0.6 


0,4 


0.2 


0.0 


50.     100.    150.    200.    250.    300 


S-MELLIN  FREQUENCY 


Shape  1  Jr   Features 
Figure  20 
88 


a. 


200. 


S-MELLIN  FREQUENCY 


300 


Shape  2  FM  Features 
Figure  21 
89 


effect   is   noted   with  respect   to  the  difference  squared 

analysis   in   Table  2t.  Although   TRI2   is   rore   closely 

correlated   to  TRIl  than  RECT2  the  squared  error  shows  just 

the  reverse.  The  results  are  encouraging.  The   preprocessor 

has  greatly   simplified  the  classification  problem  for  the 
canonic  shapes  above. 

2.   Ship  Signatures  and  Results 

A  single  ship  was  used  to  make  these  preliminary 
tests  for  this  thesis.  Signatures  were  taken  every  ten 
degrees  around  a  ship  from  zero  to  fifty  degrees.  The 
results  are  plotted  and  compared  over  twenty  degree  aspects 
in  figures  22-23.  The  signatures  are  the  result  of  very 
high  resolution  radar  signature  data  that  has  been  degraded 
and  smoothed  to  a  lower  resolution  with  essentially  no 
noise  present.  Recalling  that  the  purpose  of  the 
preprocessor  was  to  remove  variance  due  to  pure  shifting 
and  pure  scaling,  leaving  enough  information  for 
classification,  to  the  eye  there  seems  to  be  little 
encouragement  from  these  results.  It  is  recalled  that  a 
goal  of  this  preprocessor  is  to  make  the  classifiers  job 
easier  by  removing  dependence  on  shifting  and  scaling  of 
the  original  data.  The  information  that  remains  depends  on 
unspecified  signal  characteristics  that  here  appear  to 
useful  in  discriminating  shape  classes   and   possibly   ship 
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classes.  However  no  real  conclusion  can  be  drawn  at  this 
point  because  of  the  sirall  data  base  and  the  absence  of  an 
automatic  classifier  to  generate  an  optirral  feature  space. 
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Table  2 

Canonic  Shape  Fourier  -  Mellin 
Feature  Comparisons 


a.  Peak  Correlation  Values 


RECT 

RECT/2 

TRI 

TRI/2 

RECT 

1.00 

0.95 

0.50 

0.42 

RECT/2 

- 

1.00 

0.52 

0.41 

TRI 

- 

- 

1.00 

0.98 

TRI/2 

~. 

— 

— 

1.00 

b.  Squared  difference  between  features. 


RECT 

RECT/2 

TRI 

TRI/2 

RECT 

.000 

.011 

.010 

.011 

RECT/2 

- 

.000 

.015 

.014 

TRI 

- 

- 

.000 

.011 

TRI/2 

-i 

~ 

— 

.000 
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V.  CONCLUSIONS 

The  preprocessor  design  began  by  considering  a  generic 
classification  systerr.  A  distinction  was  drawn  between  the 
classifier  and  the  preprocessor.  The  preprocessor  is 
problerr  specific.  It  assists  in  the  classification  by 
extracting  a  set  of  features  for  the  classifier.  The 
extracted  feature  space  has  enhanced  natural  clustering  on 
the  basis  of  shape  by  removing  information  that  was 
extraneous  for  classification.  Two  useless  characteristics 
were  identified  as  shifting  and  scaling.  The  preprocessor 
was  designed  to  remove  dependence  on  these  two 
characteristics  by  using  the  invariant  properties  of  a 
Fourier  transform  followed  in  series  by  a  Mellin  transform. 
The  resulting  set  of  coefficients  are  Fourier-Mellin  (FM) 
features. 

A.  REVIEW 


In  Chapter  II  a  Mellin  transform  was  developed  using  the 
conventional  digital  processing  approach  which 
exponentially  warps  the  domain  and  then  transforms  the 
spectrum  by  an  FFT .  This  method  is  sometimes  identified  as 
a  fast  Mellin  transforrr  (FMT).  The  unbounded  behavior  of 
the  exponentially  warped  frequency  spectrum  was  shown  to 
result  in  a  function  that  could  not   be   transformed.   That 
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is,  the  warped  function  has  no  transform  because  the  Mellin 
integral  is  indeterrrinat  e  at  the  lower  limit.  Error 
correcting  techniques  cannot  compensate  for  the  effect 
because  the  error  itself  cannot  he  computed  in  general.  The 
hound  for  the  error  was  found  and  seen  to  he  the  envelope 
for  existing  error  correction  functions. 

Chapter  III  developed  some  useful  properties  of  the 
Mellin  transform  that  were  used  to  modify  the  signal  so 
that  the  pitfalls  isolated  in  Chapter  II  were  avoided.  This 
was  done  hy  modifying  the  input  to  the  Mellin  transform  to 
always  he  transformable.  To  simplify  the  implementation  and 
to  control  the  effects  of  sampling,  a  direct  Mellin 
transform  was  used  for  the  development  of  the  modifiers. 
Several  suitable  modifiers  were  determined  and  tested  with 
differently  scaled  inputs  for  which  the  closed  form 
solution  was  known.  The  direct  Mellin  algorithm  that 
produced  features  closest  to  the  closed  form  solution  was 
chosen  for  use  in  the  preprocessor.  With  the  modifications 
in  place,  the  preprocessor  was  tested  and  shown  to  produce 
features  that  were  invariant  to  shifting  and  scaling.  It 
was  also  shown  that  the  features  retained  enough 
information  to  classify  canonic  shapes. 

Chapter  IV  discussed  what  type  of  information  is 
required  for  classification.  Signal  structure  or  shape   was 
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identified  as  key  inf  orirat  ion .  A  discussion  of  wnat 
inforrration  is  required  for  classification  and  a  rr^eans  of 
keeping  the  signal  structure  intact  throughout  the 
preprocessor  was   advanced,    hut    not    empirically  verified. 

B.    FUTURE   WORK 

The  design  FM  preprocessor  does  produce  a  feature  space 
with  enhanced  clustering,  but  problerrs  rerrain  to  be 
resolved  before  the  full  potential  of  the  system  can  be 
realized.  There  are  three  extant  conditions  that  detract 
frorr  the  performance  of  the  implemented  preprocessor. 
First,  the  preprocessor  is  not  computationally  efficient. 
Second,  most  of  the  signal  structure  that  should  be  vital 
to  classification  is  obviously  lost.  And  third,  a  complete 
verification  of  performance  has  not  been  conducted.  The 
current  preprocessor  produces  an  enhanced  feature  base  for 
a  classifier,  but  attention  to  these  main  weak  points  will 
greatly  improve  the  applied  techniques. 

1.   Efficient  Processing 

A  direct  Mellin  transform  using  N  spectral  samples 
to  transform  into  M  Mellin  spectral  coefficients  requires  M 
by  N  multiplications.  If  only  a  few  coefficients  are  to  be 
used  then  the  number  of  multiplications  may  be  small.  Using 
the   FMT   requires   about  N(20+ln(N))  arithmetic  operations 
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for  the  Interpolation  and.  the  PFT.  If  twenty-five  or  more 
Mellin  coefficients  are  required,  the  FMT  is  faster.  The 
trodlfiers  used  to  prepare  the  spectrurr  for  the  direct 
r^ellin  transform  will  also  work  for  as  FMT,  but  this  should 
te  demonstrated  eiperirrentally .  Because  the  FMT  directly 
weights  the  lower  numbered  samples  it  may  produce  more 
accurate  results  as  well. 

2.   Conservation  of  Information 

Chapter  IV  discussed  the  type  of  information 
required  for  visual  pattern  recognition  in  humans.  The 
preservation  of  this  information  should  be  a  specified 
design  goal  for  the  preprocessor.  The  preprocessors  built 
for  this  thesis  removed  much  of  these  vital  signal 
properies.  A  means  was  introduced  to  limit  or  control  the 
loss  of  structural  detail  by  zeroing  the  fundamental  phase 
and  adjusting  each  of  the  remaining  complex  coefficients  to 
reconstruct  phase  relationships.  This  may  be  done  in  the 
frequency  domain  as  described,  or  by  a  similar  operation  in 
the  time  domain,  shifting  the  centroid  to  zero.  In  either 
case,  to  transform  more  information  about  the  signal  will 
effectively  Increase  the  sensitivity  of  the  features  to 
characteristics  in  the  time  domain.  This  increase  in 
sensitivity  needs  measured  to  confirm  the  approach.  It  is 
also   possible   that    continued   selective    zeroing   of 
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coefficients  may  offer  iirproved  perforrrance  or  robustness 
to  the  system  as  a  vhole. 

3 .   Verification 

Although  the  improved  performance  was  demonstrated 
with  respect  to  ealier  ft^  digital  preprocessing,  a  direct 
irrprovement  factor  needs  to  be  established.  Time  domain 
correlation  should  be  used  as  a  measure  of  original  signal 
classif lability .  This  treasure  needs  to  be  compared  to  the 
F^  dorrain  correlation  recorded  as  Table  2  in  Chapter  IV. 
Next,  the  preprocessor  or  an  improved  version  will  have  to 
be  married  to  a  classifier  and  realistic  data  used  to 
evaluate  its  effect  on  the  classification  system.  The 
preprocessor  built  was  design  to  be  used  on-line.  An 
on-line  classifier  needs  to  be  built  as  well. 

The  systematic  evaluation  of  ship  profiles  using  F^. 
features  is  still  a  requirement.  For  several  ships,  IM 
features  must  be  singled  out  and  plotted  together  as  a 
function  of  aspect  angle.  The  purposes  are  to  establish  a 
range  of  aspect  over  which  classification  may  be  possible, 
to  evaluate  changes  of  structural  content  as  discussed  in 
section  two  above,  and  to  determine  the  beam  signature  as  a 
classification  "node".  Although  these  purposes  assist  in 
the  system  design,  the  third  is  more  of  an  operational 
necessity.  All  ship  signatures  will   degrade   to   the  beam 
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aspect  node'  so  that  this  hypervolume  in  the  feature  space 
is  occupied  in  corrmon.  Therefore  the  bearr  condition  rust  he 
detected  and  withheld  from  ever  entering  the  classifier. 
Nor  should  beaii'  signatures  be  used  for  classifier  training. 
The  heam  "node"  needs  to  he  deteririned  separate  from  the 
classifier.  Initially  this  information  can  be  plotted  to 
eiarrine  feature  behavior  and  confirm  the  methodology, 
automated  methods  will  quickly  follow.  These  analysis 
functions  may  never  reside  in  the  classification  system 
itself,  but  must  be  a  part  of  the  tools  used  to  develop  a 
workable  classification  system. 
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APPENDIX  k 


C 

C  *  THIS  PHOGRAM  IS  A  DIGITAL  IMPLEMENTATION  OF  * 

C  *  A  FASI  FOOHIiia  TRANSFORM  FOLLOWED  BY  A  aELLiN  * 

C  ♦  TRANSFORM.   THE  DIGITAL  IMPLEMENTATION  OF  THE  * 

C  *  BELLIN  IS  AS  AN  EXPONENTIALLY  SAMPLED  SPECTRUM  * 

C  *  IHEN  RON  THROOGH  AN  FFT  AND  CORRECTED  FOR  THE  * 

C  ♦  ERROR  BY  ANY  ONE  OF  SEVERAL  CORRECTION  ♦ 

C  ♦  SOBROOTINES,  CORCTX.  * 

C 

C 

C         RECORDED  PLOTS  /  PLOT  N0MEE2  OSING  RECORD  CALLS 

C  -TIME  FUNCTION  (CONT)  /  1 

C  -SAMPLED  TIME  FUNCTION  /  1 

C  -FFT  (CONT)  /  2 

C  -EXPONENTIALLY  SAMPLED  FFT  /  2 

C  -UNIFORMLY  SAMPLED  FFT  /  3 

C  -DISCRETE  MELLIN  FEATURES  /  4 

C 

C 

C      RECORDED  PLOTS  /  PLOT  NUMBER  USING  RECANG  CALLS 

C  -TIME  FUNCTION  (CONT)  /  1 

C  -FFT  UNIFORMLY  SAMPLED  /  2  (MAG  &    PHASE) 

C  -MELLIN  FEATURES  /  3  (MAG  &  PHASE) 

C 

C 

C 

DIMENSION  XSEAL (200) ^XIMAG (200)  ,XINT (200) ,1  (64)  , 

ISCALE  -    31 

au=5 

M=2**Ma 

NU=7 

N=2**N0 

CALL  SAMP  (XREAL-XIMAG.N) 

CALL  RECORD  fXREAL, XIMAG , N, ISCALE) 
C 

C      COMPUTE  THE  ACTUAL  SAMPLE  POINTS  AS  PROVIDED  3Y  THE 
C      SAMPLED  VIDEO. 
C 

N0=5 

N=2**N0 


CALL  SAMP  (XEEAL.XIMAG.N) 
CALL  RECORD  (XRELl,XIMAG,N, IS 


SCALE) 


ISCALE  =  31 

LU=7 

L=2«*LU 

EC  200  1=1, L 


IF(I.GT.N)  GO  TO  150 
XINT  (I)=XREAL(I) 
PRT  (I)=XIMAjG(I) 


GO  TO  2  00 
150    XINT  (I) =0.0 

PRT  (I)=0.0 
200    CONTINUE 
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CALL    FPT(XINT,PRT,L,LO) 


CALL    PFT(XBEAI*,1IMAG,N,NU) 
310         CALL    HECOED  (XR2AL,XIilAGtN,ISCALE) 
CALL    BECOBD ixINT,PHT,L,ISCALE) 

,E  =  ;•' 

10    1= 

:no£ 


ISCALE    =31 

DO    300    1=1, N 

XBEAL  (I)=SQaT  (XBEAL  (I)  **2+XIHAG  (Ii**2) 

XIMAG  (I)=0 
300  CCNTI! 

C 
C 

C  CALCDLATE    THE    NEW    SAMPLE    POINTS     AND    INTEBPOLATE 

C  IC    FIND    THE    EXPONENTIALLY    SAMPLED    VALUES. 

C 

CALL  NOPTS  (XINT,M,N) 
C 

C      AFTES  NEW  INTJ^RVALS  CALCULATED  AND  STORED  TEMPOEABILY 
C      IN  VECTOR  XINT^  STORE  FOB  LATEB  PRINTING  IN  T  VECTOR. 


C 


DO  400  1=1, M 
PBT  (I)  =0.0 

^  (i)=xir-  " 


I  (I)=XINT  (I) 
UOO    CCNTINOB 
C 

C      HITH  THE  NEH  TIMES  IN  VECTOR  XINT,  AND  THE  CURRENT 
C      SPECTRUM  SA.MPLES  IN  VECTCB  XBEAL,  COiiPUTE  THE  NEW 
C      EXPONENTIAL  SAMPLES  AND  ENTER  THESE  INTO  VECTOR 
C      XINT  BY  USING  THE  CHOSEN  INTERPOLATION  METHOD. 
C      FINALLY-  RECORD  THE  FIRST  SAMPLE  TO  BE    USED  LATER 
C      TO  CORRECT  THE  MELLIN  TRANSFORM  FOR  LOW  FREQUENCY 
C      LOST  IN  THIS  EXPONENTIALLY  SAMPLED  TECHNIQUE. 
C 

F0=XREAL(1) 
399    CONTINUE 

CALL  INTP2(XHEAL,N,XINT,M) 

WRITE  {2,401)M 
401    FORMAT  (14) 

CALL  SMAX (XINT, PRT,M, SCALE) 

CO  410  1=1, M 

XREAL  (I)=XINT(I) 

XINT(I)  =SCALE  *  XINT  (I) 

XIMAG  (I)=0 

WRITE (2,415)TiI) ,XINT  (I) 
415    FCRMAT(?10.5,7X,F10.5) 
410    CONTINUE 
C 

C      SUBMIT  THESE  EXPONENTIALLY  SAMPLED  VALUES  TO  THE 
C      FINAL  FFT  BLOCK. 
C 

CALL  FFT (XBEAL, XIMAG, M,MU) 

C      APPLY  THE  COBRECTION  TERM,  FIND  THE  MAGNITUDE 

C      OF  THE  NOW  SCALE  AND  TIME  INVARIENT  FEATURES. 

C 

CCC         CALL    C0RCT1  (XREAL, XIMAG, M,F0) 

CALL    C0RCT2    XREAL, XIMAG, M,FOj 

CALL    RECORD  (XREAL, XIMAG, M,ISCALE) 

STOP 

END 
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c 
c 
c 
c 
c 
c 


«« :|i  «4t «  41 «  4i  4>«  *  «  Hi^' 4c  ♦«  4c  *  *  4t9<t  *  «  «  4i4(4(  ♦  3|(]k«  :ti « !4i:4c:K:tc :»:  lie « 
*    SUBROUTINE    FFT,    GIVEN    THZ   COMPLEX  « 


%&. 
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101 


100 


♦  N  C0EFFI2CIENTS  USING  A  CFT .  * 

SUBROUTINE  EFT (XHEAL, XIMAG- N ,N0) 

DIMENSION  XBEAL(N)  ,XIMAG  (N) 

N2=N/2 

N01=ND-1 

K=0 

DO  100  L=1,N0 

DO  101  I=l-N2 

E=IBIIR  (K/2**NU1.NU) 

AfiG=6.2&3185*P/FLOAT(N) 

C=COS  <ABG) 

S=SIN  (ABG) 

K1=K+1 

TEEAL=XHEAL(K1N2)  ♦C+XIMAG  (K1N2)  *S 

TIMAG=XIMAG (K1N2)  *C-XBEAI (K1N2) *S 

XBEAL  (K  1N2)  =XfiEAL  (K  1)  -TBEAL 

XIMAG  (K1N2i  =XIMAGiK1)-TIMAG 

XEEAL  Ik  1)  =XBEAL  (KT)  +TREA1 

XIMAG  (K  1)  =XIMAG  (K1)  -HIM AG 

K=K>1 

K=K>N2 

IFJK.LT.N)  GO  TO  102 

NU1=NU1-1 

N2=N2/2 

DC  103  K=1,N 

I=IBITE  (K'T.NO) +1 

IF(I.LE.K)GC  TO  103 

TREAL=XBEAL  (K] 

TIMAG-XIMAGJK 

XBEAL  (K)=XB£A] 
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C 

C 

c 
c 


200 


I 


XIMAG  (K 
XBEAL (I 
XIMAG  (l' 
CCNTINU: 
BETURN 
END 


^XIMAG 

=TR£AL 
=TIilAG 


a) 


lEITR 
IBITR(J,NU) 


FUNCTION 

FUNCTION 

J1=J 

IEITR=0 

DO  200  1=1, NO 

J2=J1/2 

I£ITH=IBITE*2+ (J1-2*J2) 

J1=J2 

BETURN 

END 
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C  *  INTP2  IS  A  SECOND  OfiDEfi  INTERPOLATION  BASED  ON  A  * 

C  *  PCLINOMIAL  EQUATING  TO  FIEST  AND  SECOND  DEfilVATIVES  * 

C  *  APPROXIMATED  BY  CENTRAL  DIFFERENCES.  THE  INPUT  VECTOR  * 

C  ♦  <X>  CONTAINS  UNIFORMLY  SAMPLED  DATA  WITH  UNITY  ♦ 

C  *  BETWEEN  CONSECUTIVE  SAMPLES.   VECTOR  <Y>  IS  INPUT  * 

C  *  WITH  THE  NEW  SAMPLE  TIMES  AND  OUTPUT  WITH  * 

C  *  THE  NEW  SAMPLES.  THERE  ARE  N  UNIFORM  SAMPLES  IN  <X>  * 

C  *  AND  M  WARPED  SAMPLES  IN  <Y>.  * 

C 

c 

SUBROUTINE    INTP2 (X. N, Y, M) 

DIMENSION  X3(3)  ,X(N)  ,Y(M) 
C 
C 

C      CHOSE  THE  X  SAMPLE  TIME  CLOSES  TO  THE  WARPED 
C      TIME  HELD  IN  <Y> 
C 

DO  60  IY=1,M 

ET  =  100 

DO  10  IX=1.N 

IXTIME  =  IX  -  1 

CTABS  =  ABS  (DT) 

CIST  =  FLOAT  (IXTIME)  -  Y  (lY) 

DIST  -    ABS(DIST) 

IF  (DTABS  .LT.  DIST)  GO  TO  25 

DI  =  Y(IY)  -  IXTIME 
10     CONTINUE 
25     IXTIME  =  IXTIME  -  1 

J1  =  -1 

DO  30  J=1,3 

J1  =  J  •••  IXTIME  -  1 

X3(J)  =  X  (Jl) 
30     CGNTINOE 

TIME  =  Y(IY) 

Xai)    '    FCN^(X3,  TIME, IXTIME) 
60     CONTINUE 

RETURN 

END 
C 

C      FCN  IS  THE  INTERPOLATION  RULE. 
C 

FUNCTION  FCN(X,T,I) 

DIMENSION  113) 

TX  =  FLOAT  (I) 
C 

C      COMPUTE  THE  COEFFICIENTS. 
C 

A  =  (X(3)  -  2.*  X(2)  *    X(1))  /  2. 

E  =  JXJ3)  -  Xri))  /  2.0  -  2.*  A  *  TX 

C  =  X(2)  -  A  *  IX**2.  -  B  *  TX 
C 

C  CCMPUIE    THE    INTERPOLATED    VALUE. 

C 

FCN    =    A    ♦     (T**2)     *    B   *    1    ■¥    C 

RETURN 

END 
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C  *  SOBfiOOTINE    INTEfiP    USES    A    LAGRANGE    THIfiD  * 

C  *  CBDEB    METHOE    Ofi    A    SECOND    OEDEH    POLYNOaiAL  * 

C  *  IC    COMPUTE    IHE    INPUT    SAMPLE    HAVEFORM.  * 

C  ♦  * 

C  ♦  THE    SUBROUTINE    INPUTS    ARE:  * 

C  *  X-SAMPLED    HAVEfOEM  ♦ 

C  *  N-THE    NUMBER    OF    X    SAMPLES  * 

C  *  Y-INPUT    SAMPLE    TIMES  ♦ 

C  *  -OUTPUT    ITERPOLATED    SAMPLES  * 

C  *  M-THE    NUMBER    OF    I    SAMPLES  ♦ 

C 

c 

SUBROUTINE  INTERP  (X  ,N,  Y  ,  MJ 
EIMENSION  Xa(4)  ,X(N)  ,Y  (fl)  ,SHFTX  (156) 

C 

C  FIRST    CHOOSE    XTIME    NEAREST    EACH    YTIME, 

C  THEN    COMPUTE   THE    INTERPOLATED   VALUE    AT    THAT    PT. 

C 

c 

CO    100    1=1, N 
15=1+5 


IF(I.GE.  (N-5)  )     GO    TO    101 
-  I5)=X(I) 

GO    TO    100 


SHFTX{I5] 


101         15=1- (N-5) 

SHFTx!l5)  =X(I) 
100    CONTINUE 

CO  40  1=1, M 
UO     CONTINUE 

EO  60  IY=1,a 

YP=100 

EC  10  IX=1,.N 

IXT=IX-1 

YPABS=AES  (YP) 

DIST=ABS(IX.I-Y(IY)  ) 

IF(YPABS.LT.DIST)  GO  TO  25 

YP=Y(IY)-IXI 

IXU=IX+5 
10  CONTINUE 
25     IX4=IXU-2 

DO  30  J=1,4 

X4{J)=SHFTX(IX4+J) 
30     CONTINUE 

IM1=IX4>1 

10=1X4+2 

IEl=IX4+3 

IP2=IX4+4 
35     Y  (IY)=YLAGR(X4,YP) 
60     CONTINOE 

RETURN 

END 


ins 


c 

C      FONCTION  YLIN  MAKES  k   LINEAR  INTERPOLATION. 

C 

C 

FUNCTION  YLIN(X4,yP) 

ClflENSION  X.U(4) 

IF(YP.LT.O.)  GO  TO  50 

YLIN=XU  (2)  ♦YP*  (XU  (3)-X4  (2) ) 

BETORN 
50     YIIN=X4  (2)-YP*(X4(1)-X4  (2)) 

RETURN 

£N0 
C 

c 
c 

C  fONCTION  YLAGH  COMPOTES  THE  LAGRANGE  MULTI-PLIEES 

C  AND  MAKES  THE  INTERPOLATION  FOR  A  CHOSEN  OFFSET 

C  FROM  THE  CENTRAL  SAMPLE  X4(2) 
C 

c 

FUNCTION  YLAGR(X4,YP) 
DIMENSION  X4(4) 
CM1=-YP*iYP-T)*(YP-2)/6.0 
C0=  (YP**2-1)* (YP-2) /2.0 
CP1=-YP*  (YP  +  1)  *  (YP-2)  /2-Q 
CP2=YE*  (YP**2-1)/6 .  0 

YLAGR=CM1«X4(1)  4-C0*X4  (2)  +CP1*X4  (3)  ♦CP2*X4  (4) 
BETORN 
END 
C 

C  ♦  SOBROOTINE  NOPTS  CALCULATES  THE  M  EXPONENTIAL  SAMPLE  * 

C  ♦  POINTS  FROM  THE  N  UNIFORM  SAMPLES  OF  THE  EXISTING     ♦ 

C  ♦  SPECTROM  IN  PREPARATION  FOR  AN  INTERPOLATION.   THIS   ♦ 

C  ♦  EXPONENTIALLY  SAMPLED  SET  OF  POINTS  ARE  STORED  IN     ♦ 

C  *  THE  INPOT  VECTOR  X.  * 

C 


SOBROOTINE  NOPTS (X,K,N) 

DIMENSION  X.(156) 

0N=PLCAT(N)/2.0  +  1.0 

£M=FLCAT(M)-1.0 

DELZ=LOG(0N)./EM 

DC  100  1=1, 1j 

SI=FLOAI(I)-1.0 

VALUE  =  SI*DELZ 

X(IV  =  EXP  (VALUE) 


100    CONTINOE 
RETURN 
END 
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C         *    THIS    IS    A    STDB    PROVIDING    k    TEST    SKIN    BETOBN  * 

C         ♦    TC    SiaOLATE    A    SKIN    RETUfiN.       NORMALLY    THIS    SUBROUTINE    * 
C         *    aiLL    ACTUALLY    SAiJPLE    A    SKIN    fiETORN.  * 

C 

SUBROUTINE  SAMP  (XR, XI, N) 
CIHENSICN  XB(N)  ,XI  (N) 

C      DESIGNATE  SCALE  FACTOR  AND  UNSCALED  TIME  SHIFT. 


C 


SCALE  =  16.  /  32. 

SHIFT  =  O.C  /  32. 

SCALE  =  1 .0/SCALE 

T0=  .5  ■•■  saiFT 

N1=N-1 

EQ  100  1=1, N 

XR(I)=0. 

21 


m°o: 


C 

C      CALCULATE  THE  TIME  OF  THE  SAMPLE. 
C 

SK=  (FLOAT  (I) -1.0) /FLOAT  (N1) 
C 

C      EUILC  THE  TEST  SKIN  RETURN. 
C      IN  THIS  (C  2-PER)  CASE  A  DOUBLE  PEfilMID. 

TSCALE  =  (SK-TO)*SCALE 

IHFM  =  1 

IF(IHFM  . 

IF(IWEM  . 
10     H  =  8.  / 

IF  (TSCALE 

IF  (TSCALE 


IF  JTSCALE 

[I)  =  i 

GO  TO  100 


XE(I)  =  fTO 


20  XR(I)  =  H  ♦  10.0  -  (TSCALE  -  TO)  *  10.0 

GO  TO  100 

50  CONTINUE 

C  IHS  FOLLOWING  HAVE  FORM  IS  A  SQUARE  HAVE  16  SAMPLES 

C  HIDE  CENTERED  AROUND  THE  SIXTEENTH  SAMPLE. 

C  SCALING  AND  SHIFTING  DONE  ABOVE  HILL  EFFECT  THE 

C  HAVEFCHM  ACCORDINGLY. 

TSCALE    =     (SK-TO)*SCALS 
EDGE    =    1.    /    4. 

IF     (TSCALE    .LT.    -EDGE)     GO   TO    100 
IF     (TSCALE    .GT.    +EDGE)     GO   TO    100 
1B(1)    -    1.0 
100         CCNTINUE 
RETURN 
END 
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C  *  IHIS  SOBfiCUTINE  RECORDS  A  SELECTED  DATA  SET,  * 

C  ♦  AND  SCALES  THE  INDEX  TO  32  SAMPLES  (0-31).  * 

C  *  THE  INPUTS  ARE  XREAL.  XIMAG,  AND  N  JTHE  * 

C  *  NOMBEB  OF  SAMPLES) .  ISCALE  IS  A  CCNTEOL  * 

C  ♦  VARIAELE.  * 

C  *       IF  ISCALE  =  0  NO  SCALING  IS  PERFORMED  * 

C  ♦               SCALE  =  1,  AXIS  =31  * 

C  *               SHIFT  =  1  (STARTS  AT  T  =  0  )  * 

C  ♦          ISCALE  <  0  SHIFTING  AND  SCALING  OCCURS  * 

C  *                IISCALEJ  =  AXIS  LEN.  A  SHIFT  IS  INCORP-  ♦ 

C  *               ERATED  SO  THAT  AXIS  STARTS  AT  1 .  ♦ 

C  ♦          ISCALE  >  0  AXIS  SCALING  OCCORS  BUT  THERE  IS   * 

C  *               NO  SHIFT  (IE.,  THE  AXIS  STAETS  AT  0.)    * 

C  ♦  * 

C 

SUBROUTINE  RECORD (XE, XI -N, I SCALE) 

DIMENSION  XB(N)  ,XI(N)  ,RIc(200) 

SCALE=1.0 

AXIS  =  31 

SHIFT  =  1 

WRITE  (2,50)  N 
50     FORMAT (la) 

IFJISCALE  .EQ.  0)  GO  TO  40 

AXIS  =  FLOAT (ISCALE) 

AXIS  =  ABS(AXIS) 

IF(  ISCALE  .LT.  0)  SHIFT  =  0 

CALL  SMAX rXR, XI, N, SCALE) 
40     DO  100  I=T,N 

SI=  (I-SHIFT)*AXIS/ (FLOAT  (N)  -1.0) 

RECJI)=SQRT(iI(I)**2+XR  (I)**2) 
60     CONTINUE    \ 

REC  (I) =SCALE*REC (1) 

WRITE  (2,7  5)  SI -RECri) 
75     FOEMAI(F10.5,7X,FT0.5) 
100    CONTINUE 

RETURN 

END 
C 

c 


SUBROUTINE    SMAX (XR, XI, N , SCALE) 
DIMENSION    XE(N)  ,XI  (N)  ,T  (200) 
XMAX=1.0 
DO    100    I=T.N 

T  (I),=SCfiTiXR(I)  **2  +  XI  (I)**2) 
IF(T(I)  .LT.XMAX)     GO    TO    100 
XMAX    =    Tri) 
SCALE=1.0/XMAX 
100         CONTINUE 
RETURN 
END 


108 


C  *  IHIS  CORRECTION  SUBROUTINE  USES  ONE  OF  * 

C  ♦  THE  SIMPLER  CORRECTIONS  FOR  THE  MELLIN  * 

C  *  TRANSFORM.  THE  CORRECTION  IS  A  PURE  IMAGINARY.   * 

C  ♦  * 

C  *     CORRECTION  =  -FO/OMEGA  * 

C  *  ♦ 

C  *  THEN  A  MODIFICATION  IS  MADE  TO  THE  ENTIRE  * 

C  ♦  TRANSFORM  BY  A  MULTIPLICATION  BY  OMEGA.  * 

SUBROUTINE  C0RCT1 (XR, XI ,N,FO) 
CIMENSICN  XS(N;  ,XI(N} 

DO  100  1=1, N 

OMEGA  =  FLOAT (I)  -  1. 

XR(I)  =  XR(I)  *  OMEGA 
IF  (I  .EQ.  1).G0  TO  90 

XI  (I)  =  (XI  (I)  -  FO/OMEGA)  ♦  OMEGA 
GC  TO  100 
90     XI  (I)  =  FO 

100    CONTINUE 

RETURN 
END 


C  *  THIS    CORRECTION    APPLIES    THE    MORE    COMPLECATED  * 

C  *  EXPRESSION,  * 

C  ♦  CORRECTION    =    FO/2    ■••    JCOT  (FO/OMEGA)  * 

C  ♦  * 

C  *  AND    THEN    MODIFIES    THE    ENTIRE    TRANSFORM    3Y    1/OMEGA.  * 

C 

SUBROUTINE  C0RCT2 (XR, XI ,N,FO) 
DIMENSION  XR(N),XI(N) 

EC  100  1-1, N 

CMEGA  =  FLQAT(I)  -  1. 

XR{I)  =  (XR(I)  *  FO/2.)  *  OMEGA 

IF(I  .EQ.  1)G0  TO  90 

XI  (I)  =  (XI(I)  -  (FO/2.  )*COTAN  (OMEGA)  )  ♦OMEGA 

GO  TO  100 

90     XI(I)  =  FO/2. 

100    CONTINUE 

RETURN 
IND 
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APPENDIX  B 


C 

C  *  THIS  PROGEAM-  FOURIER  DIRECT  MELLIN,  TAKES  AN  INPUT    * 

C  ♦  WAVEFORM  FROM  LOGICAL  DEVICE  2,     PEEFORaS  AN  FFT        * 

C  *  FOLLOWED  BY  A  flOMT  CUTPUTING  THE  FEATURES  TO  LOGICAL   « 

C  *  DEVICE  3  FOR  LATER  PLOTTING  BY  MELPLT.   THE  MAJOR      * 

C  ♦  DATA  STRUCTURES  AND  SUBROUTINES  ARE  LISTED  BELOH       * 

C  *  IN  ORDER  OF  THEIR  APPEARANCE.   THE  SUBROUTINES  ARE     * 

C  *  DESCRIBED  IN  MORE  DETAIL  WHERE  THEY  ARE  ACTUALLY       * 

C  *  LISTED  IN  THE  PROGRAM-                                   * 

C 
C 

c 

c 

C  HAJOB  DATA  STRUCTURES: 

C  <HFM>  -  THE  INPUT  WAVEFORM  (REAL) 

C  <XFM>  -  THE  MELLIN  TRANSFORM  (COMPLSi) 

C  <CPHI>  -  REAL  MELLIN  COEFFICIENTS 

C  <SPHI>  -  IMAGINARY  MELLIN  COEFFICIENTS 

C  <STAND>  -    AN  ARRAY  HELD  FOR  LATER  COMPARISON 

C  OR  OTHER  USE. 

C  <PRI>  -  AN  ARRAY  NORMALLY  USED  TO  HOLD  REAL 

C  DATA  TEMPORILY.  A  WORK  SPACE, 

C  <IXT>  -  X  AXIS  TITLE  FOR  PLOTTING 

C  <IYT>  -  Y  AXIS  TITLE  FOR  PLOTTING 

C  <KEY>  -  NUMBER  OF  PLOT  THIS  GRAPH 

C 

c 
c 

C  SUBROUTINES: 

C  WAVE  -  READS  AN  ARRAY  FROM  LOGICAL  DEVICE  2, 

C  AND  FILLS  ZEROS  TO  MAKE  A  TOTAL  OF  256 

C  SAMPLES.  THE  OUTPUT  IS  IN  <WFM>. 

C  FFT  -  AH  FFT  BLOCK 

C  COEF  -  COMPUTES  THE  MELLIN  TRANSFORM  SAMPLE 

C  WEIGHTS.  THESE  ARE  COMPLEX 

C  NUMBERS  WHOSE  REAL  AND  IMAGINARY 

C  PARTS  ARE  STORED  IK  <CPHI>  AND 

C  <S?HI>  RESPECTIVELY. 

C  DMTM  -  APPLIES  A  MODIFIED  DIRECT  MELLIN 

C  TRANSFORM  TO  AN  INPUT  WAVEFORM 

C  PUTTING  THE  OUTPUT  IN  <XFM>. 

C  THE  ALGORITHM  IS  BASED  ON  A  FIRST 

C  BACKWARD  DIFFERENCE. 

C  SMT  -  APPLIES  A  MODIFIED  MELLIN  TRANSFORM 

C  BASED  ON  A  SECOND  DIFFERENCE. 

C  SMT2  -  APPLIES  A  MODIFIED  MELLIN  TRANSFORM- 

C  DIFFEBENT  THAN  SMT.  BUT  ALSO  BASED 

C  ON  THE  SECOND  DIFFERENCE. 

C  CDMT  -  APPLIES  A  MODIFIED  AELLIU    TRANSFORM 

C  JUST  AS  DMTM  ABOVE,  EXCEPT  THAI  THE 

C  CENTRAL  DIFFERENCE  IS  USED. 
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C  Siaf  -  APPLIES  A  MODIFIED  MELLIN  TRANSFOHN 

C  USIBG  A  BACKHASD  DIFFEfiZNCE  AS  IN 

C  DHTM,  EXCEPT  THAT  THE  INTEGRATION 

C  IS  BY  SIMPSON'S  RULE  INSTEAD  OF  THE 

C  TRAPEZOIDAL  ROLE. 

C  XAB  -  TAKES  THE  MAGNITODE  OF  THE  COMPLEX 

C  TRANSFORM  <XFM>  AND  POTS  THE 

C  MANITGDE  IN  A  SPECIFIED  VECTOR. 

C  STOK  -  NORMALIZES  A  VECTOR  BY  ITS 

C  MAGNITUDE  AND  WRITES  II  TO 

C  LOGICAL  DEVICE  3  WITH  A  TITLE 

C  FROM  LOGICAL  DEVICE  i*. 

C  HOLD  -  LOADS  ONE  VECTOR  INTO  ANOTHER. 

C  ALTER  -  CHANGES  <aFM>  BI  SCALE  S/OR  SHIFT 

C  AND  OOTPOTS  TO  A  SPECIFIED  ARRAY 

C  INTP3  -  A  SECOND  ORDER  SPLINE  INTERPOLATION, 

C  CFOHM  -  PROVIDES  TWO  CLOSED  FORM  SOLUTIONS 

C  FOR  VERIFYING  THE  MELLIN  ALGORITHMS. 

C  TITLE  -  ENTITLES  THE  PLOTS  ON  THE  BASIS  OF 

C  THE  CALLING  PROGRAM. 

C 

c 
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SC  THE  MAIN  PROGRAil  STARTS! 

DIMENSION  P:BT(256)  ,  WFM  (256)  , STAND  (256)  ,lf  (5) 
COMMON  XFM  (256.2)  -CfHI  (256, 12b)  ,SPHI (256, 128) ,PI, 

iixT(iO)-,iYT(i©f  ,Ki;y 

DATA  IF/«   PR«,'2Q0£'r' NCI  •,•     »,»     •/ 
PI  =  3.  141592654 


C      HOH  MANY  WAVEFORMS  ARE  TO  BE  TRANSFaRMED? 

BEAD(2,  10: 
10     FORMAT  (14^ 


BEAD  (2, 10)  NUMHFM 
) 


C      NTMS  IS  THE  NUMBER  OF  TIME  SAMPLES  (INCLUDING 

C      ANY  ZERO  FILLING) .   IT  IS  A  POWER  OF  TWO 

C      FOR  THE  CONVENIENCE  OF  THE  FFT. 

C      HPTS  IS  THE  NUMBER  OF  SAMPLES  INPUT  TO  THE 

C      MELLIN  TRANSFORM  BLOCK.   THE  COEFFICIENTS 

C      ABE  COMPUTED  NOW. 

NU  =  8 

NTMS  =  2**N'U 

MPTS  =  NTMS/2 
50     FCRMAI(I4)^ 

NCOEF  =  NTMS 

CALL  CO EF  (NCOEF, MPTS) 

C      SET  UP  THE  LOOP  FOR  THE  NUMBER  OF  WAVEFORMS 
C      TO  BE  PROCESSED. 

DO  500  IWAVE=1,NUM»FM 

C      GET  THE  NEXT  INPUT  WAVEFORM. 
CALL  WAVE  (WFM, NTMS) 

C      ZERO  THE  <STAND>  VECTOR  TO  BE  USED  AS  THE 
C      IMAGINARY  PART  OF  THE  NTMS  TIME  SAMPLES. 

DO  100  I=1,HTMS 

STAND  (I>  =  Q.O 
100    CONTINUE 

CALL  STOW fWFM, NTMS) 

C      TAKE  THE  FFT. 

CALL  TITLE(IF) 

CALL  FFT (WFM,STAND, NTMS, NC) 

C      TAKE  THE  MAGNITUDE  AND  PUT  IT  INTO  THE  COMMON 
C      WAVEFORM  <WFM>. 

DO  200  1=1, NTMS 

WFM(I)  =  WFM(I)**2  +  STAND(I)**2 

WFMjl)  =  SCBT(WFM(I)) 
200    CONTINUE 

CALL  STOW  (WFM,NTMS) 

C  TAKE    THE    MEiLLIN    TRANSFORM   OF    THIS    SPECTRUM 

C  USING    THE    FIRST    HALE   OF    THE    FFT    SAMPLES, 

C  OTHERWISE    KNOWN    AS    MPTS=NTMS/2.        THESE    ARE    THE 

C  ONLY    ONIQUE    VALUES. 

CALL    DMTM (WFM, MPTS, NCOEF) 

CALL    XAE(PRT,NCOEFf 

CALL    STOW  (PBT, NCOEF) 
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c 
c 

C  NEXT  CALL  THE  SECONC  OBDZE  RULE  DMT 

C  SUBROOTIilES.   BOTH  COMPOTE  THE  MELLIN 

C  USING  THE  SECOND  DIFFERENCE  APPROXIMATION 

C  INSTEAD  OF  THE  FIRST  DIFFERENCE  APPROXIMATION 

C  AEOVE.   OTHERHISE  THE  APPROACH  IS  THE  SAME, 

C 

c 

CALL    SMI(HFM,MPTS^NCOEF) 
CALL    XAEiPRT.NCOEFJ 
CALL    STCH  (PfiT,NCOEF) 

CALL    SMT2  (HPM,MPTS,NCOEF) 
CALL    XAB(PB.T,NC0EF1 
CALL    STCi3(PBT,NC0EF) 

C  CALL    CDMT    H-HICH    USES   THE   CENTRAL    DIIFERENCE 

C  ROLE    FOR    APPROXIMATING    THE    TRANSFORM. 


CALL    CDMT  (HFM,MPTS,NC0EF) 
CALL    XAB 
CALL    STCJ 


CALL    XABfPRT.NCOEF) 
;»  (PRT,NCOEF) 


C  CALL    SIMP    WHICH    COMPUTES    THE    MELLIN    TRANSFORM 

C  USING    THE    FIRST    DIFFERENCE    ALGORITHM    AND 

C  SIMPSON'S    ROLE    TO    COMPUTE   THE    MODIFIED 

C  MELLIN    TRANSFORM. 

C 

CALL  SIMP  (WFM,MPTS,NCOEF) 
CALL  XAB(PRT,NCOEF) 
CALL  STCW  (PST,NCOEP) 

C  END  THE  TRANSFORM  LOOP.  THE  TRANSFORM  HAS  BEEN 

C  OUTPOT  TO  LOGICAL  DEVICE  3  AND  PREPARED  MIIH 

C  TITLE  INFORMATION  PROVIDED  BI  LOGICAL  DEVICE  2 

C  FOE  PLOTTING  IIITH  PROGRAM  MELPLT  FORTRAN. 

C  STAY  IN  THE  LOOP  II    MORE  iiAVEFORMS  ARE  AVAILABLE. 

500  CONTINUE 

CALL    CFOEM (PRTrNCOEF) 
CALL    CF0RM(PRT,NCO£F) 

STOP 
END 
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c 
c 

C  ♦    IHE    «AVE    SOBROaTINE    PRODUCSS    A    HAVEPOfiW-    IN  * 

C  ♦    BEAD    PROM    TBE    LOGICAL    UNIT    2,    NOfifiALLY    A    SHIP  * 

C  ♦    DATA    FILE.       THE    N    SAMPLES    ARE    READ    AND    THEN  * 

C  ♦    ZEROS    ARE    STUFFED    TO    Fill   THE    256    SAiiPLES.  ♦ 

C  ♦    THE    REQUESTED    BAVEFORil    IS   OUTPUT    IN    THE  * 

C  ♦    COMMON    ARRAY    HFM  (256) .  * 

C 

SUBROUTINE  BAVE (HFM ,NPTS) 
DIMENSION  WFfl(N?TS) ,ID(5[ 

COMMON  XFM  (256,2)  ,CPHI  (256, 128)  ,SPHI  (256, 128) ,PI, 
IIXT(IO)  ,IYT(10)  ,KEy 
DATA  ID/'  RAN«,«GE  S« , • AMPL • , ' ES   *,«     •/ 

READ (2, 10) KEY 

READ  (2,  30)  (IXI  (I)  ,1  =  1,10) 

READ(2,30)  IYT(I>  ,1=1,10) 
30     FORMAT ( 10 A^) 
50     HEAD (2,  10) N 

READ(2,20i.  (B?M(I)  ,1  =  1, ») 

CALL  TITLE  (ID> 
10     FORMAT  (14) 
20     FORMATJFIO.5)^ 

IF(N  .EQ.  NPTS)  RETURN 

N  =  N  >  1 

DO  100  I=N,NPTS 

«FM(I)  =  0.0 
100    CONTINUE 

BETURN 

END 
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c 
c 
c 
c 
c 
c 
c 
c 
c 


*  SUBHOOTINE    PFT,    GIVEiJ    THE   COMPLEX  * 

*  SAMPLES     (2**Na=N) ,     WILL    HETUEN       THE  * 


USI 


♦    H    COEfFIECIENTS    USING    A    £FT .  * 


102 


101 


100 


103 


C 
C 

c 
c 
c 
c 


200 


SUBBCDTIN 
DIMENSION 
N2=N/2 
N01=ND-1 
K=0 

DO  ICO  L=1,NU 
CO  101  I=1'n2 
E=IBITRiK 
ABG=6.283 
C=COS  (ABG 
S=SIN(ABG 
1=K+1 


E    EFT(XHEAI-XIMAG,N,NO) 
Xfi£AL(N)  ,XIMAG  (N) 


i:n2 

/2**NU1,NU) 
18'5*P/FL0AT(N) 


K 
K1N2=K1*N 

tseal=xbe 
timag=xim 
xbeal  (kin 
ximag (kin 

XBEAL  (K1) 
XIMAG  (K1) 
K  =  K-»'1 
K=K+N2 

ifVk.lt.n 

N01=N01-1 

N2=N2/2 

DO  103  K= 

I=IBITE  (K 

IF/I.LE.K 

TBEAL=XBE 

TIMAG=XIM 

XBEAL 

XIMAG 

XBEAL 

XIHAG 

CCNTINO! 

BETUBN 

END 


AL  (K1N2)*C+XIMAG 
AG  iKlN2)*C-XBEAI 
2)=XBEAL(K1)-TREAL 
2)  =XIMAGiK1j -TIMAG 
=XBEAL(K1) ♦TREAL 
=  XIMAG(K1)  -►TIMAG 


)  GO  TO  10  2 


I.N 

-1-NU) -H 
)GC  TO  103 

Alr(K) 

AGJK) 
XBEAL  (I) 
XIMAG (I) 

TBEAL 
TIMAG 


(K1N2)  *S 
iK1-N2)  *S 


IBITB 
IBITE(J,NU) 


FUNCTION 

FUNCTION 

J1=J 

IBITB=0 

DO    200    1=1, NU 

J2=J1/2 

IBITR=IBITB*2^-  (J1-2«J2) 

J1=J2 

BETOBN 

END 
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c 

C  *  THE  COEF  SUBROUTINE  COMPUTES  THE  ilELLIN  * 

C  ♦  COEFFICIENTS  IN  TO  COMMON  ARBAYS  CPHI  AND  * 

C  *  SPHI  THAT  REPRESENT  THE  REAL  AND  IMAGINAiiY  * 

C  ♦  FARTS  RESPECTIVELr.  THE  TERMS  ARE  COMPUTED  * 

C  ♦  EY  THE  FORMULA:  * 

C  ♦     PHI  (I, J)  =  J**S  ,  »H£RE  S  =  A  NORMALIZED  * 

C  *  DISCRETE  RADIAN  FREQUENCY.  * 

C 

SUBROUTINE  COEFfNCOEF, NPTS) 

COMMON  XPM(256,2)  ,CPHI(256,  128)  ,SPHI  (256,  128)  ,PI, 
lIXTilO)  -lYTMO)  ,KEY 

CO  TOO  I  =  1,NC0Ef 

RI  =  FLOATill 

CMSGA  =  2,0   *  PI  *  RI  /  36. 

DO  20C  J  =  1,NPTS 

RJ  =  FLOAT  (J) 

CPHI  (I,  J)  =  COS  (OMEGA  *  ALOG  (RJ)  ) 

SPHI  (I- J)  =  SIN  (OMEGA  *  ALOG  (RJ)  ) 
200    CONTINUE 
100    CONTINUE 

BETURN 

END 
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c 
c 

Q        «4i«:(i]|i4t4i4i4i4i4i;|iik4(«*4(«««**4t«4t*4i]|c*4(«*3|t*i((*«i|ij(e**«4(:it;|i:4tiK:4i«4c4i«4c4i4c 

C  *  IHE  DMTM  SOBfiOUTINE  PEEFORMS  A  DESCRETE  MELLIN  ♦ 

C  *  TRANSrORil  ON  THE  ASfiAI  HfM.   THE  FORMULA  * 

C  ♦  FOR  ONE  MELLIN  FREQUENCY  VALUE  IS:  * 

C  *       XFM  (I)=SUM(K=1  10  NETS)  (HFM  (1+1) -WFil(I)  )  *K**S  * 

C  ♦  * 

C  *  THE  COMPLEX  COMPONENTS  OF  K**S  ARE  COMPUTED  PRIOR  * 

C  *  TO  CALLING  DMTM  AND  STORED  IN  THE  COMMON  AHiiAXS  * 

C  ♦  CPHI  AND  SPHI  AS  REAL  AND  IMAGINARY  PARTS  * 

C  *  filSPECTIVELY.  THE  AIGOEIIHM  IS  BASED  ON  THE  * 

C  ♦  FIRST  DIFFERENCE.  THE  TRAPEZOIDAL  RULE  IS  USED  * 

C  *  FOR  THE  INTEGRATION.   THE  COMPLEX  OUTPUT  FOR  * 

C  ♦  THE  TRANSFORM  IN  IN  THE  COMMON  ARRAY  <XFM>.  ♦ 

C 
C 

SUBROUTINE  DMTM (SAMP, NPTS^NCOEF) 
DIMENSION  SAMP(NPTS) ,ID (5) 

COMMON  XPM  (256,2)  ,CPHI(^56, 128)  ,SPHI  (256, 128) ,PI, 
1IXT  (10)  ,IYI  (10)  .KEY 
DATA  ID/«D-ME»,'LLIN»,«  FHE • , * QUEN» , * CY   •/ 

CALL  TITLE  (ID) 
DC  100  I=1,NC0EF 
XFM  (I,1)=0.0 
XFM  a-2j=0.0 
SI  ==  NPTS  -  1 
CO  200  J=1,N1 
J1  =  J  +  1 

XFM  (1,1)  =  CPHI  (I,  J)  *  (SAMP(J>  -  SAMP(J1))  ♦  XFM  (1,1) 

XFmJi,2)_  =  SPHlil,J)  *  iSAMP(J)  -   SAMP(JI))  >  XFM  (1,2) 

200    CON-flNUE        »    /    *  \    »    t 

100    CONTINUE 
RETURN 
FND 
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c 

C  ♦  SMT  IS  A  SUBROUTINE  THAT  COMPUTES  A  NUMERICAL  * 

C  *  APPROXIMATICN  TO  THE  MEILIN  TRANSFORM  AS  DOES  * 

C  *  EMTM  ABOVE.   SMT  USES  A  TRAPAZOIDAL  APPROXiaATION  * 

C  ♦TO  COMPUTE  THE  TRANSFORM,  BUT  USES  THE  SAME  * 

C  *  COEFFICIENT  MATRIXES  <CPHI>  AND  <S?HI>  CONTAINED  ♦ 

C  *  IN  COMMON.  » 

c 


SUBROUTINE    SMT  (SAilP  ^NPTS  ,NCOEF) 
DIMENSION    S.AMP(NPTSi,IDi5) 

COMMON    XFM  (256,2)  ,CPHI  (z56, 128)  , SPHI (256,12  8) ,PI, 
1IXT(10I  ^lYT  (10)  .KEY 
DATA    ID/» S-ME», 'LLIN«,»     FRE« , 'QUEN* , 'CI       •/ 
CALL    IITLE(ID) 

C  INITIALIZE    IHE      INPUT    ARRAY    AND    COMPUTE 

C  THE    LOOP    CONSTANTS. 

N1    =    NPTS    -    1 

N2    =    N1    -    1 

C  SET  UP  THE  TRANSFORM  LOOP.  THE  OUTER  LOOP 

C  SETS  UP  THE  COEFFIECIENTS  WHILE  THE  INNER 

C  LOOP  COMPUTES  THE  SOM  WHICH  ARE  THE 

C  COEFFICIENTS. 


^x«.o,  ,)  =  Q.O 
XFM  (J,2i  =  Q.O 


DO  200  J  '    1,NC0EF 
XFM(J,1 

XFM Jj.r 

DC  TOO 

10  =  I 

11  =  I  ♦  1 

12  =  I  ♦  2 

DELTA  =  SAMP(IO)  -  2.*  SAMP  (II)  +  SAMP  (12) 
XFa(J,1)  =  XFM<J,1)  ■«•  DEITA*CPHI  (J,I)*I 


XFa(J,1)  =  XFM(J,1)  ■«•  DEITA*CPHI  (J,I)*I 
XFMJJ-2i  =  XFM  (J, 2)  ♦  DELTA*SPHI  (J,I) *I 
100    CONTINUE 

C      XFM(J,1)  =  XFM(J,1)  /  (FLOAT  (J)  *PI/18.) 
C      XFM  (J, 2)  =  XFM  (J,  2)  /  (FLOAT  (J)  *PI/18. ) 

XFM(J,1)  =  XFH(J,1)  /  SQRT(1*(FL0AT(J)*PI/18.)**2) 
XFM(J,2)_  =  XFM  (J, 2)  /  SQRT  ( 1+ (FLOAT  (J)  *PI/1  8.)  **2) 
200    CONTINUE 
BETURN 
END 


118 


c 

c 

C  *  SMT2  IS  A  SOBSOOTINE  THAT  COMPaTES  A  NUMEfilCAL  ♦ 

C  *  AFPSOXIMATICN  TO  THZ  MEillN  TRANSFOEil  AS  DOES  * 

C  *  EMTM  ABOVE.   SMT  USES  A  TEAPAZOIDAL  APPSOXIiUATION  * 

C  *  TO  COMPUTE  THE  TRANSFOBM-  BUT  USES  THE  SAME  * 

C  *  COEFFICIENT  MATRIXES  <CPHI>  AND  <SPHI>  CO^JTAINED  ♦ 

C  *  IN  CO:iflON.  * 

C 

SDBSOOTINE  SMT2 (SAMP, NPTS,NCCEF) 
DIMENSION  SAMP (NPTS) .10(5) 

COMMON  XFM  (256,2) ,CPHI (256, 128)  , SPHI (256,128)  ,PI, 
IIXT(IO)  ,IYT  (10)  .KEY 
DATA  ID/»S2-M«, 'SLLI' ,• N  FH • , *EQUE» ,^ NC Y  •/ 
CALL  TITLE  (ID) 

C      INITIALIZE  THE   INPUT  A25AI  AND  COMPUTE 
C      THE  LCOP  CONSTANTS. 

N1  =  NPTS  -  1 

N2  =  N1  -  1 

C  SET  UP  THE  TRANSFORM  LOCP.  THE  OUTER  LOOP 

C  SETS  UP  THE  COEFFIECIENTS  WHILE  THE  INNER 

C  LCOP  COMPUTES  THE  SUM  HHICH  ARE  THE 

C  COEFFICIENTS. 

DC  200  J  =  1,NC0SF 
XFM(J,1)  =  0.0 
XFM  JJ,2)  =  Q.O 
EO  100  I  =  .1,N2 

10  =  I 

11  =  I  -«•  1 

12  =  I  +  2 

DELTA  =  I*  (SAMP  (10)  -  2.*  SAMP  (in  ♦  SAMP  (12)) 
EELTA  =  DELTA  +  (SAMP  (12)  -  SAMP  (10)_) 
XFM(J,1)  =  XFM(J,1)  *    DELTA^CPHI  (J,I) 
XFMjJ-2i  =  XFM  (J, 2)  *    DELTA*SPHI  (J,I) 
CONTINOE 


100 


200 


REMOVE    THE    COLOfi 


XFM(J,1)  =  XFH(J,1)  /  SQfiT(U(FLOAT  (J)*PI/18.)**2) 
XFMJJ,2j  =  XFM  (J, 2)  /  SCET  ( 1  +  (FLOAT  (J)  *PI/18. )  **2) 
CONTINUE 


RETURN 
END 
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c 
c 

C  *  THE  CDMT  SUBROUTINE  FERfCEMS  A  DESCHETE  MELLIN  » 

C  *  TRANSFORM  ON  THE  COMMON  ARRAY  WFM.   THE  FORMULA  * 

C  *  FOR  ONE  MELLIN  FREQUENCY  VALUE  IS:  * 

C  *f  XFM  (I)=SnM(K=1  TO  NPTS)  (MFM  (I*  1) -«FM  (I- 1)  )  *K**S  * 

C  *  THE  COMPLEX  COMPONENTS  OF  K**S  ARE  COMPUTED  PRIOR  * 

C  *  TO  CALLING  CDMT  AND  STORED  IN  THE  COMMON  ARRAYS  * 

C  *  CPHI  AND  SPHI  AS  REAL  AND  IMAGINARY  PARTS  * 

C  *  RESPECTIVELY.  THE  CENTRAL  DIFFERENCE  IS  USED.  * 

c 
c 


SUBROUTINE  CDMT (H, NETS , NCOEF) 
DIMENSION  HJNPTS) ,ID(5) 

COMMON  XFM  (256,2)  , CPHI  (256, 128)  ,SPHI  (256, 128)  ,PI, 
1IXT(10)_,IYT(10)  -.KEY 
DATA  ID/»C-aE«  r'^I'LIN',  »  FRE  •  ,  *  QUEN*  ,  •  CY   •/ 
CALL  TITLE  (ID) 


C      SET  UP  DERIVATIVE  OF  INPUT  VECTOR  <H (I)  > 
C      IDENTIFIED  HERE  AS  <G>. 

NG  =  NPTS  -  1 

DO  200  J  =  1  ,  NCOEF 
XFM  (J,1)  -  Q.O 
XFM(J,2)  =  0,0 

OMEGA  =  FLCII(J)  *  PI  /  18. 

EC  100  I  =  1,NG 

IM1  =  I 

IP1  =  I  +  2 

Dfl  :=  (H  (IP1)  -  H(IM1))  /  2. 

C      COMPUTE  THE  J-TH  COEFFICIENTS  BY  THE  SUM. 

XFM(J,1)  =  XFM(J,1)  +  DH  ♦  CPHI  (J, I) 
XFM  (J, 2)  =  XFM  (J, 2)  ♦  DH  *  SPHI  (J, I) 

100    CCNTINOE 
C 

COLOR  =  1 

XFM(J,1)  =  XFM(J,1)  ♦  COLOR 

XFM  (J, 2)  =  XFM  (J,  2)  *  COLOR 

200    CONTINUE 

BETURN 
IND 
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C  *  THIS    SOBflOOTINS    USES    SIMPSON'S    RULE    TO    APPROXIMATE  * 

C  ♦  THE    MODIFIED    MELLIN    TRANSFORM.       THE    MODIFICATION  * 

C  «  IS    FREQUENCY    TIMES    THE    EREQUENCY    DERIVATIVE    OF  * 

C  ♦  THE    FFT.  * 

C 


SDBBOUTINE  SIMP (H, NETS, NCOEF) 
CIMENSION  HiNPIS) ,ID(5) 
1  XFM  {256,2)  .CPHI 

DATA  II);"I-riE«',  »LLIN«,  •  ERE  •  , 'QUEN*  ,  *  CY   •/ 


COMMON  XFM  (256,2)  ,CPHI  (256, 128)  ,SPHI (256, 128)  ,PI, 
IIXT(IO)  ,IYT(10)  -KEY 


CALL  TITLE  (ID) 


N1  =  NPTS  -1 
N2  =  NPTS  -  2 


DO  100  J=1,NC0EF 

FSIMP  =  2.0 

CMEGA  =  PI  *  FLOAT  (J)  /  18, 

COLOR  =1 


XFM  (J,1)  =  0. 
XFM  (J, 2)  -    0. 


DC  200  1=1, N1 

IM1  =  I 

10  =  1  +  I 

IF  (FSIMP  .GT.  3.)  GO-  TC  67 

FSIMP  =  4. 

GO  TO  69 

67 

ESIMP  =  2. 

69 

CONTINUE 

DELIA  =  H(I0)  -  HjIMI) 
XFM(J,1   =  XFM(J,T   ♦  FSIMP 
XFM  (J, 2)  =  XFM  (J, 2)  >  FSII^P 

♦  DELTA  ♦ 

*  DELTA  ♦ 

200 

CONTINUE 

100 

CONTINUE 

RETURN 

END 

CPHI(J,IM1) 
SPHI(J,IM1) 
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c 

C  *  THE  XAB  SOEBOUTINE  TkKES    THE  MAGNITUDE  OF  THE  * 

C  «  COMMON  COaFLEX  ABEAX  XFt5  (256,2)  AND  PLACES  * 

C  ♦  THESE  VALUES  IN  OUTPUT  VECTOfi  <XMAG>  FOB  LATEfi  * 

C  «  EBINTING.  ♦ 

c 
c 

SUBROUTINE    XAB (XHAG ,NPTS) 
DIMENSION    XflAG(NPTS) 

CCaaON    XFM  (256. 2) .CPHI (256, 128)  ,SPHI (556, 128) ,PI, 
IIXTMO)  ,IYT(10f  ,kSY 
00    too    I=1,NPTS 

XMAGri)    =    SQRT(XFM  (I,  1)  ♦♦2+XFM  (I,2r**2) 
100         CONTINUE 
BETUBN 
END 
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c 

C  *  THE  STOW  SOBEOUTINE  STORES  THE  INPUT  AfifiAI  Pfil(NPTS)  ♦ 

C  ♦  INTO  THE  LOGICAL  DEVICE  2  FOB  LATER  USE  IN  PLOTTING.  * 

C  ♦  NPLOl   NUMBEfiS  THE  PLOTS  FOE  LATER  IDENTIFICETION,  * 

C  *  AND  A  PLOT  TITLE  FOR  THE  tfELLIN  FREQUENCY  IS  ADDED  ♦ 

C  *  FOR  CONVENIENCE.  ♦ 

C  *        THE  MODULUS  OF  THE  * 

C  *  TRANSFORM  TAKEN,  SCALED  TO  UNIT  MAGNITUDE,  AND  ♦ 

C  *  OUTPUT  TO  LOGICAL  UNIT  2.  * 

C  ♦     INPUT:     PRI  -  TO  BE  SCALED  TO  1  AND  WRITTEN  » 

C  ♦  TO  LOGICAL  UNIT  2.  * 

C  *  NPLOT  -  THE  NUMBER  OF  THE  PLOT  * 

C  ♦  KEY  -  NUMBER  OF  CURVE  THIS  PLOT  * 

c 


SUBROUTINE  STOH (PfiT,NPTS) 
DIMENSION  P5T(NPTS) 

COMMON    XFMJ256,2)  , CPHI (256, 128)  , SPHI (256,128)  ,PI, 
IIXT(IO)  ,IIT(1G)  /key 

WRITE  (3,13) NPTS 

WRITE  (3,1  3)_ KEY 
CKK         IF     (KEY    .NE.    ^lGO    TO    12 
C  WRITE    AXIS   XABELS 

WRITE  (3, 10)  (IXI  (I), 1=1,  10) 

WRITEJ3.IO)    lYT(I)  *  1=1,10) 

10         forma1(1oa4)  ' 

12        nelot  =     npiot  +  1 

fiPET  =  0.0 
13     FORMAT  (14) 

DO  100  I=1,NPTS 

IF  (PET  (I)  .GT.  BPRT)  BFET  =  PET  (I) 
100    CCNTINOE 

IF  (BPRT  .LI.  .00001)  BPET  =  1. 

DO  200  I=1,NPTS 

PET  (I)  =  PRT(I)  /  BPET 

I  =  FIOATil) 

WRITE  (3-20)T,PRT(I) 
20     FOaMAT(2(3X,f9.5)) 
200    CCNTINUE 

RETURN 

END 
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c 
c 

C  *  THE    HOLD    SOBfiOOTINE    TAKES    THE    INPUT    FILE  * 

C  ♦  <FILIN>    AND    STORES    IT    IN    THE    OUTPUT    FILE  * 

C  *  <fILOUT>    FOB    TEMPORARY    STORAGE.     TfiE    FILE  * 

C  ♦  <FILIN>    REMAINS    UNCHANGEE.  ♦ 

C 

SUBROUTINE  HOLD (FILIN,FILOUT,NPTS) 
CIMENSION  FILIN (NPTS) ,FILOUT  (NPTS) 

DO  10C  I  =  1,NPTS 
FILOUT(I)  =  FILIN  (I) 
100    CONTINUE 
RETURN 
END 

C 

C  *  THE  ALTER  SUBROUTINE  »IIL  ALTER  THE  COiiMON  ARRAY  * 

C  ♦  <HFM>  AS  SPECIFIED  BY  THE  INPUT  VARIABLES  * 

C  *  <SCALE>  AND  <SHIFT>.   THE  ALTERED  WFM  IS  OUTPUT  * 

C  *  IN  THE  VECTOR  <ALT>.  * 

c 
c 

SUBROUTINE  ALTER (ALI, »FM , SCALE, SHIFT, NPTS) 
DIMENSION  ALTiNPTS)  .«FM  (NPTS) .TOLD  (256)  -TNEH  (256) 
COMMON  XFM  (256,2)  ,C§HI (^56, 128)  ,SPHI  (256, 128)  ,PI, 
lIXTilC)  ,IYt  (10)1  , KEY 
DO  loo  1=1, NPTS 
TOLD  (I)  =  FLOAT  (I) 

TNEH(I)  =  TOLD  (I)  /  SCALE  +  SHIFT 
100    CONTINUE 

CALL  INTP3 (HFM, TOLD, ALT, TNEH, NPTS) 

RETURN 

END 
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c 

C  *  INTP3  IS  A  SECOND  ORDEfi  INTERPOLATION  BASED  OH  A  * 

C  ♦  ECLI2I0MIAL  EQUATING  TO  FIflST  AND  SECOND  DERIVATIVES  * 

C  *  APPROXIMATED  BI    CENTRAL  DIFFERENCES.   THE  INPUT  * 

C  *  VECTOR  <X0>  HAS  OLD  SAMPLES  AT  TIMES  IN  <I0>.  THE  * 

C  *  NEW  SAMPLE  TIMES  ARE  INPUT  THROUGH  ARRAY  <TN>  AND  * 

C  *  THE  COMPUTED  SAMPLES  AT  THESE  TIMES  ARE  OUTPUT  IN  * 

C  ♦  ARE  OUTPUT  IN  THE  VECTOR  <XN>  AND  <HFM>.  * 

C 
C 

SUBROUTINE  INTP3 (XO.TO. XN.TN.NPTS) 

TN(NPTS) 

1IXT  (10)  ,IYt  (10)  ,XEY    ■    '  '    '    •'   ' 

C 

C  CHOSE    THE    <T0>    SAMPLE    TIME    CLOSES    TO    THE    SAfiPED 

C  TIME    HELD    IN   <TN> 


C 


DC    60    I=1,NPTS 
DT   =    100 


XPTS  =  FLOAT(NPTS) 
"~((TN(I)  •GE.  1.) 
XN(I)  =  0.0 


IF((TN(I)^  .GE.  1.)  .AND.  (IN  (I)  ,LE.  XPTS))  GO  TO  5 


GO  TO  60 
5      EO  10  J=1,NPTS 

DTABS  =  ABS  (DT) 

CIST  =  TO  (J)  -  TN (I) 

DIST  =  ABS(DIST) 

IF  (DTABS  .LT.  DIST)  GO  TO  25 

DI  =  IN  (I)  -  TO  (J) 
10     CCNTINOE 
25     JTIME  =  J  -  1 

J1  =  -1 

DO  30  J=1,3 

J1  =  J  +  JTIME- 1 

IF  UJI  .GE.  1)  .AND.  (J1  .LE.  NPTS)  )  GO  TO  29 

X3(J)  =  0.0  ' 

GO  TO  30 

29  13  iJ)    -    XO  (J1) 

30  CCNTINOE 
TIMN  =  TN(I), 
IIMO  =  TO  Jj.TIME) 

XN(I)  =  FCN (X3,TIM0,TIMN) 
60     CONTINUE 

CO  500  1=1, NPTS 
500    CONTINUE 

RETURN 

END 
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c 
c 

C      PCN  IS  THE  INTERPOLATION  ROLE. 
C 

FONCTION  FCNiX,TO,TN) 

DIMENSION  X(i) 
C 

C      COMPOTE  THE  COEFFICIENTS. 
C 

A  =  1113)    -   2.*  X(2)  ^-  X(1I)  /  2. 

E  =  H\3)     -  XJ1))  /  2.0  -  2.*  A  *  TO 

C  =  X(2)  -    h   *    10**2.  -  B  *  TO 

C      COMPOTE  THE  INTERPOLATED  VALOE. 
C 

FCN  =  A  *  (IN**2)  ♦  3  *  TN  ♦  C 

BETORN 

END 
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c 

C  *  THIS  SUBROOTINE  IS  THE  EESDLT  OF  A  CLOSED  ♦ 

C  *  PORa  CALCULATION  OF  THE  CCNTINOOUS  SINC**2  ♦ 

C  ♦  EEING  TAKEN  FROM  THE  TIME  DOMAIN  THROUGH  THE  * 

C  ♦  ENTIRE  FOOBIER-ilELLIN  PROCESSING.   THE  ALGORITHM  * 

C  *  IS  USED  TO  VERIFY  THE  F2M  PROCESS.   THE  OUTPUT  ♦ 

C  *  FEATURE  SPACE  SHOULD  BE  IDENTICLE  TO  THAI  PRODUCED  * 

C  ♦BY  ANY  OF  THE  MELLIN  ALGORITHMS.   FOR  THIS  REASON  * 

C  *  THE  SAMPLE  POINTS  ARE  SYNCHRONIZED  WITH  THOSE  USED  ♦ 

C  *  ABOVE.   <NCOEF>  FM  COEFFICIENTS  ARE  USED.   THE  * 

C  ♦  COMMON  VARIABLE  <KEY>  MUST  BE  SET  TO  99  PRIOR  TO  * 

C  *  CALLING  TO  GET  THE  SINC**2  OUTPUT  FEATURE  SPACE.  * 

C  *  TO  OUTPUT  THE  CLOSED  FORM  SOLUTION  FOR  A  RAMP  IN  * 

C  *  THE  FREQUENCY  DOMAIN,  <KEY>  MUST  EQUAL  100.  * 

C 

c 

SOBBOUTINE  CFORM (CF ,NCOEF) 

DIMENSION  CE(NCOEF) 

COMMON  XFM  (256,2) ,C?HI (256, 128)  , SPHI  (256,128)  ,PI, 
IIXT(IO)  ,IYT  (10)  ,KEY 

IF(KEY  .EQ.  100)  GO  TO  200 

IF(KEY  .NE.99)  RETURN 

DO  100  M=1,NC0EF 

XM  =  FLOAT  (M)  -  1.0 

HMAG  =  1  +  (PI  ♦  XM  /  18.)**2 

CF(M)  =  SQRT(HMAG)  /  HMAG 

100  CONTINUE 
GO  TO  101 

200    DO  202  I  =  1,NC0EF 

XI  =  FLOAT  (I)  -1.0 

OMEGA  =  PI  ♦  XI  /  18. 

CF(I)  -    2./SQRT(4.  *   0MEGA**2) 
202    CONTINUE 

101  CONTINUE 
BEAD (2, 20) KEY 


READ(2^10)  (IXT  (I)  ,1  =  1,10) 
READ  (2,  10)  lYT(I)  ,1  =  1,10) 

10     FORMAT (10A4) 

20     FORMAT  (I^) 


CALL  STOW  (CF,NCOEF) 

RETURN 
END 
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c 

C  *  TITLE  TITLES  THE  I  AXIS  FOB  THE  PLOTTING  * 
C  ♦  EBOGBAa  ACCORDING  TO  THE  ALGORITHa  aSED  TO  * 
C    ♦  GENEBATE  THE  fEATOBE  SPACE.  ♦ 

SOBHOUTINE  TITLE  (ID) 
DIMENSION  IDiS) 

COMMON  XPM  (256,2) ,C2HI  (256, 128)  ,SPHI (256, 128) ,PI^ 
lIXTilO)  ,IYTi10)  ,KEY  v    #    /#   # 

DO  too  1=1, D 
lYT(I)  =  ID(I) 
100    CONTINOE 
EETUBN 
END 
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